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ABSTRACT 


Below  a  specific  threshold  signal-to-noise  ratio  (SNR),  the  mean  squared  error  (MSE)  perfor¬ 
mance  of  signal  direction-of-arrival  (DOA)  estimates  derived  from  the  Capon  algorithm  degrades 
swiftly.  Prediction  of  this  threshold  SNR  point  is  of  practical  significance  for  robust  system  design 
and  analysis.  The  exact  pairwise  error  probabilities  for  the  Capon  (and  Bartlett)  algorithm  are 
derived  herein,  given  by  simple  finite  sums  involving  no  numerical  integration,  include  finite  sample 
effects,  and  hold  for  an  arbitrary  colored  data  covariance.  An  accurate  large  sample  approximation 
of  these  error  probabilities  in  terms  of  the  well  tabulated  complementary  error  function  is  also 
provided.  Via  an  adaptation  of  an  interval  error-based  method,  these  error  probabilities,  along 
with  the  local  error  MSE  predictions  of  Vaidyanathan  and  Buckley,  facilitate  accurate  prediction 
of  the  Capon  threshold  region  DOA  MSE  performance  for  an  arbitrary  number  of  well  separated 
sources,  circumventing  the  need  for  numerous  Monte  Carlo  simulations.  A  large  sample  closed  form 
approximation  for  the  Capon  (and  Bartlett)  threshold  SNR  is  provided  for  uniform  linear  arrays. 
A  new  exact  two-point  measure  of  the  Capon  probability  of  resolution,  that  includes  the  deleterious 
effects  of  signal  model  mismatch,  is  a  serendipitous  by-product  of  this  analysis  that  predicts  the 
SNRs  required  for  closely  spaced  sources  to  be  mutually  resolvable  by  the  Capon  algorithm.  Lastly, 
a  new  general  strategy  is  provided  for  obtaining  accurate  MSE  predictions  that  account  for  signal 
model  mismatch. 
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1.  INTRODUCTION 


The  threshold  region  mean  squared  error  (MSE)  performance  of  angle  estimates  derived  from 
the  Capon  algorithm,  a.k.a  the  minimum  variance  distortionless  response  (MVDR)  spectral  esti¬ 
mator,  is  the  primary  subject  of  this  report.  The  classical  Bartlett  algorithm  wherein  an  angle 
estimate  is  obtained  via  a  smoothed  periodogram  is  a  closely  related  approach  whose  analysis  shall 
be  developed  in  parallel.  The  Capon  algorithm  can  be  described  as  a  high-resolution ,  data  adaptive, 
nonparametric  approach  to  direction  of  arrival  (DOA)  estimation  [15,  17,  18,  33,  29,  11,  52,  61], 
whereas  the  Bartlett  algorithm  is  a  classical,  nonparametric,  conventional  beamforming  approach 
[6,  7,  14,  29,  52].  The  parallel  development  given  herein  is  motivated  by  the  similarities  involved 
in  the  analysis  as  well  as  the  contrast  provided  by  use  of  conventional  beamforming  performance 
as  a  reference  point  for  super-resolution  approaches.  The  angle  estimates  are  obtained  as  the 
arguments  of  the  largest  peaks  of  an  estimated  power  spectrum  over  spatial  frequency  or  angle. 
Similar  to  maximum-likelihood  (ML)  estimation,  the  Capon  and  Bartlett  algorithms  are  beamscan 
type  algorithms  involving  a  nonlinear  search  of  a  data  dependent  objective  search  function  (OSF). 
Parameter  estimates  derived  from  nonlinear  searches  often  exhibit  a  threshold  effect  in  MSE  perfor¬ 
mance  [66,  59,  43,  65].  Below  a  specific  signal- to- noise  ratio  (SNR)  called  the  estimation  threshold 
point,  the  MSE  departs  from  an  often  well-behaved  (inversely  proportional  to  SNR)  asymptotic 
MSE  performance  and  rises  rapidly.  It  rises  until  it  reaches  a  maximum  that  at  times  can  be  well 
approximated  by  the  worst  case  variance,  Le.,  that  of  an  estimate  that  is  uniformly  distributed  over 
the  angle  search  domain.  The  SNR  at  which  the  MSE  performance  achieves  this  level  of  futility 
is  called  the  no  information  point.  Figure  1  illustrates  this  composite  MSE  performance  typical  of 
nonlinear  estimation  schemes.  Three  distinct  regions  of  MSE  performance  are  evident:  no  infor¬ 
mation  (as  SNR— ►  0),  threshold  (or  ambiguity),  and  asymptotic  (as  SNR— ►  oo).  The  departure  of 
the  MSE  curve  from  the  asymptotic  performance  as  the  SNR  decreases  below  the  threshold  point 
is  due  to  global  errors  caused  by  the  presence  of  false  peaks  (sidelobes)  in  the  underlying  ambiguity 
function1.  Each  scan  over  the  angle  search  domain  of  interest  of  the  OSF  produces  a  realization 
of  a  stochastic  process  whose  statistical  behavior  ultimately  determines  the  accuracy  of  the  an¬ 
gle  estimates.  The  ambiguity  function  has  a  strong  influence  on  the  character  and  nature  of  this 
non-Gaussian,  nonstationary  stochastic  process,  especially  in  the  vicinity  of  the  threshold  region. 
The  composite  MSE  behavior  shown  in  Figure  1  is  typical  of  nonlinear  ML  estimation  [59,  60],  but 
likewise  occurs  with  the  Capon  and  Bartlett  DOA  estimators.  Although  the  threshold  effect  of  the 
Capon  algorithm  has  been  observed  in  practice  [61],  accurate  prediction  of  the  composite  MSE  per¬ 
formance  curve  is  an  open  problem.  Because  the  Capon  algorithm  is  a  high-resolution  approach, 
significant  beamsplit  ratios  (BSR)  are  possible  in  the  threshold  region,  making  its  performance 
prediction  in  this  region  of  practical  interest.  The  goal  of  this  analysis  is  to  provide  accurate  pre¬ 
diction  of  the  threshold  SNR  point  for  angle  estimates  derived  from  the  Capon  algorithm,  as  well 
as  prediction  of  the  threshold  region  of  the  MSE  curve,  accounting  for  finite  sample  effects  [16,  35], 
and  colored  noise.  Signal  model  mismatch  is  known  to  be  a  formidable  practicality  to  realizing 

lrFhe  ambiguity  function  for  the  Capon  estimator  is  obtained  by  evaluating  the  OSF  over  the  angle  range  of 
interest,  but  using  the  true  data  covariance  in  place  of  the  usual  estimated  data  covariance.  Although  an  ambiguity 
function  is  classically  defined  as  the  deterministic  signal  component  of  the  output  of  a  linear  matched  filter  response 
[(56,  46,  44  ,  60],  it  is  appropriate  to  define  the  ambiguity  function  of  the  Capon  algorithm  as  stated,  given  the  stochastic 
nature  of  the  signal  models  considered  herein. 


1 


Figure  1.  Typical  composite  MSE  curve  for  nonlinear  parameter  estimation. 

the  full  potential  of  super-resolution  approaches  [18,  26,  63,  64].  Thus,  some  consideration  is  given 
herein  to  account  for  its  presence  in  these  threshold  region  MSE  predictions. 

A  classical  method  of  MSE  approximation,  referred  to  herein  as  the  method  of  interval  errors 
(MIE),  was  introduced  by  Van  Trees  [59]  and  provides  a  means  of  predicting  threshold  region  per¬ 
formance  of  nonlinear  estimation  techniques.  MIE  (and  its  variants)  is  a  well  established  approach 
that  has  been  shown  to  provide  accurate  MSE  prediction  of  ML  estimation  and  related  nonlinear 
techniques  well  into  the  estimation  threshold  region  (see  discussion  in  Section  2).  Although  MIE  is 
well  established,  obtaining  the  components  necessary  to  apply  this  approach  to  a  specific  algorithm 
can  be  quite  nontrivial.  Indeed,  the  successful  application  of  MIE  requires  good  approximations  of 
two  quantities:  (i)  interval  error  probabilities,  and  (ii)  the  asymptotic  MSE  performance.  Both  of 
these  quantities  are  algorithm  dependent.  The  interval  error  probabilities  quantify  the  likelihood 
that  the  estimator  derives  its  signal  parameter  estimate  from  an  interval  (local  neighborhood)  of  the 
search  domain  dominated  by  a  false  peak  of  the  ambiguity  function,  as  opposed  to  the  local  interval 
containing  the  true  peak  (global  maximum).  These  probabilities  are  often  well  approximated  via 
the  Union  Bound  in  conjunction  with  good  approximations  of  pairwise  error  probabilities.  Exact 
pairwise  error  probabilities  for  the  Capon  (and  Bartlett)  estimator  are  derived  herein  that  hold  for 
an  arbitrary  colored  data  covariance  and  include  finite  sample  support  effects.  The  Cramer-Rao 
Bound  (CRB)  often  serves  as  an  adequate  approximation  of  the  asymptotic  local  error  MSE  perfor- 
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mance  of  ML  estimation  due  to  its  well  established  properties  regarding  efficiency  [59].  Although 
the  Capon  estimator  is  not  asymptotically  efficient  [61],  its  asymptotic  (in  SNR  and  large  sample 
support)  local  error  MSE  performance  has  been  studied  extensively.  Stoica  et  al.  [51]  established 
the  theoretical  large  sample  MSE  performance  of  the  Capon  estimator  as  well  as  an  improved 
Capon  estimator  in  the  context  of  DOA  estimation.  Vaidyanathan  and  Buckley  (VB)  [57,  58]  in¬ 
dependently  provided  theoretical  predictions  of  this  MSE  performance  that  includes  the  effects  of 
finite  training  and  sensitivity  to  random  perturbations  of  the  assumed  signal  model.  Hawkes  and 
Nehorai  (HN)  [21]  have  extended  these  results  to  the  estimation  of  vector  signal  parameters  and 
include  a  parallel  theoretical  development  for  the  Bartlett  spectral  estimator.  The  MSE  predictions 
of  VB  that  capture  finite  sample  effects  shall  be  used  herein  for  the  Capon  algorithm,  and  those  of 
HN  shall  be  used  for  the  Bartlett  algorithm. 

Signal  modeling  errors  are  inevitable  in  several  application  areas.  Assessing  the  impact  of 
such  errors  on  the  accuracy  of  DOA  estimates  is  of  practical  concern.  A  new  strategy  is  proposed 
for  applying  MIE  to  cases  involving  signal  model  match.  It  consists  of  a  modification  of  this 
classic  technique.  The  proposed  strategy  is  based  primarily  on  the  fundamental  nature  of  nonlinear 
searches  driven  by  multimodal  ambiguity  functions,  and  it  represents  a  unique  generalization  of 
this  classic  approach. 

A  serendipitous  by-product  of  this  present  threshold  region  analysis  is  the  definition  of  a  new 
exact  two-point  measure  of  the  Capon  probability  of  resolution  that  provides  accurate  prediction  of 
the  SNRs  required  for  closely  spaced  sources  to  be  mutually  resolvable.  It  is  obtained  via  a  modified 
pairwise  error  probability  calculation,  holds  for  an  arbitrary  colored  data  covariance,  and  includes 
finite  sample  effects.  This  resolution  measure  is  parameterized  by  the  true  data  covariance  and 
choice  of  scanning  vectors.  It,  therefore,  can  account  for  the  presence  of  signal  model  mismatch, 
allowing  for  sensitivity  analyses  of  signal  modeling  errors  [20.  64].  The  approach  taken  herein  leads 
to  an  exact  probability  density  function  (pdf)  for  the  ratio  of  two  Capon  (and  Bartlett)  power 
spectral  estimates,  from  which  the  error  probabilities  easily  follow. 

This  present  analysis  provides  two  very  useful  results  in  relation  to  the  asymptotic  local  error 
MSE  approximations.  First,  for  signals  that  are  spaced  by  more  than  a  beamwidth,  this  analysis 
provides  accurate  prediction  of  the  threshold  SNR  point,  thus  indicating  the  region  of  SNRs  for 
which  these  asymptotic  MSE  approximations  are  valid.  Second,  for  closely  spaced  sources,  the 
proposed  two-point  measure  of  the  Capon  algorithm  probability  of  resolution  lower  bounds  the 
threshold  SNR  [25]  above  which  the  asymptotic  MSE  approximations  are  valid. 

This  technical  report  is  organized  as  follows.  Section  2  summarizes  the  contributions  of  early 
pioneers  of  threshold  region  MSE  analysis,  and  reviews  recent  activity  in  this  area  of  research. 
Section  3  gives  brief  descriptions  of  the  Capon  and  Bartlett  algorithms,  summarizes  known  results 
on  their  large  sample  local  error  MSE  performance  prediction,  and  establishes  the  notational  con¬ 
vention.  Section  4  describes  the  MIE  technique  of  MSE  prediction,  its  adaptation  to  the  Capon  and 
Bartlett  algorithms  for  well  separated  sources,  and  it  includes  a  stepwise  summary  for  the  computa¬ 
tion  of  the  necessary  pairwise  error  probabilities.  Section  5  proposes  a  method  of  direct  closed  form 
approximation  for  the  Capon  (and  Bartlett)  algorithm  large  sample  threshold  SNR  for  the  canoni¬ 
cal  and  theoretically  fundamental  case  of  a  single  planwave  signal  in  white  noise.  Section  6  proposes 
a  two- point  measure  of  the  probability  of  resolution  for  the  Capon  (and  Bartlett)  algorithm  that 


can  be  applied  to  the  case  of  closely  spaced  sources.  Section  7  considers  a  direct  calculation  of 
the  Capon  algorithm  resolution  SNR  (sometimes  called  the  detection  SNR  [25])  required  to  resolve 
sources  with  a  specified  level  of  confidence.  Section  8  describes  a  new  generalization  of  MIE  that 
encompasses  the  possibility  of  signal  model  mismatch.  Section  9  provides  several  numerical  results 
corroborating  the  success  and  demonstrating  the  utility  of  these  predictions.  Conclusions  are  given 
in  Section  10.  Appendices  are  provided  detailing  the  derivation  of  the  required  pairwise  error  prob¬ 
abilities  and  exact  distributions  of  the  ratio  of  power  estimates.  MATLAB  code  is  also  provided 
for  the  threshold  region  MSE  prediction  of  the  canonical  case  of  a  single  planewave  signal  in  white 
noise  for  the  interested  reader. 
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2.  THE  THRESHOLD  EFFECT:  A  HISTORICAL  PERSPECTIVE  AND 

REVIEW  OF  RELATED  WORK 


The  threshold  effect  is  in  essence  a  rapid  degradation  in  system  performance  observed  when 
the  signal  level  falls  below  a  specific  signal- tonoise  ratio  (SNR)  referred  to  as  the  threshold.  His¬ 
torically  such  an  effect  has  been  observed  for  over  half  a  century  in  the  performance  of  communica¬ 
tion  systems  that  employ  data  modulation  schemes  such  as  pulse  code  modulation  (PCM)  [45]  or 
phase/ frequency  modulation  (FM)  [32].  These  schemes,  coupled  with  a  system  designer’s  desire  for 
efficiency  that  inevitably  leads  to  a  minimum  power  requirement,  brings  threshold  region  system 
performance  to  the  forefront  of  communication  systems  analysis.  Transmission  of  binary  digits 
(called  bits)  obtained  from  a  coding  strategy  that  converts  useful  information  (speech/ voice,  for 
example)  into  a  form  robust  for  efficient  data  transmission  over  a  noisy  channel  is  the  staple  of 
all  communication  systems.  These  bitstreams  are  often  decomposed  into  sub-bitstreams  that  are 
mapped  into  one  of  several  symbols  of  a  constellation  of  say  A I  possible  symbols.  The  job  of  the 
receiver  is  to  determine  which  symbol  was  transmitted  over  the  channel  in  the  presence  of  noise.  A 
communication  system,  therefore,  has  as  one  of  its  most  fundamental  tasks  that  of  M- ary  detection, 
z.e.,  the  testing  of  multiple  hypotheses.  Correct  detection  of  each  transmitted  symbol  obviously 
permits  the  correct  synthesis  of  the  original  bitstream,  and  ultimately  intelligible  reconstruction  of 
the  original  transmitted  information.  When  the  SNR  falls  below  the  threshold,  incorrect  symbol 
detection  errors  can  corrupt  the  bitstream  making  the  received  demodulated  waveform  and  message 
(speech,  for  example)  unintelligible.  Shannon  referred  to  this  threshold  SNR  as  the  threshold  of 
intelligibility  [45]  and  offered  a  very  interesting  explanation  for  its  presence.  His  explanation  derives 
from  a  dimension  theoretical  interpretation  of  coding  theory  that  suggests  that  folds  (or  ambigui¬ 
ties)  in  the  signal  space  are  an  inevitable  topological  consequence  of  mapping  (or  coding)  a  higher 
dimensional  signal  into  a  lower  dimensional  space  or  vice  versa  when  the  low-dimensional  signal  is 
designed  to  fill  the  high  dimensional  space.  Shannon  attributes  the  threshold  effect  of  communica¬ 
tion  systems  to  the  presence  of  these  folds /ambiguities,  much  in  the  same  way  that  the  observed 
threshold  effect  in  the  MSE  performance  of  nonlinear  parameter  estimation  schemes  is  attributed 
to  the  presence  of  subsidiary  multiple  maxima  in  the  underlying  ambiguity  function  associated  with 
the  OSF. 

Woodward  [66]  draws  this  parallel  between  communication  theory  and  nonlinear  parameter 
estimation.  He  considers  the  classic  range/delay  estimation  problem  encountered  in  Radar  applica¬ 
tions  and  describes  in  lucid  detail  the  performance  accuracy  of  such  estimates.  Woodward  defines 
the  ambiguity  function  in  the  classical  sense  as  the  signal  component  of  the  optimal  matched  filter 
output  and  describes  its  impact  in  the  presence  of  noise  on  the  accuracy  of  range  estimates  for  the 
high,  moderate,  and  low  SNR  cases.  His  analysis  identifies  three  distinct  regions  of  performance 
as  SNR  is  varied2,  and  he  offers  an  interesting  interpretation  of  the  threshold  effects  in  range 
estimation  from  the  perspective  of  Shannon’s  thesis  (see  p.  90  of  [66]);  namely,  the  transmitted 
waveform  that  carries  the  information  is  a  multidimensional  quantity,  whereas  the  information  it 
carries  (the  target  range)  is  one-dimensional.  Thus,  one  can  viably  interpret  the  waveform  as  a 
multidimensional  encoding  of  the  original  message  (range)  that  is  one-dimensional.  Similarly,  one 

2  Although  Woodward  did  not  name  the  regions,  they  are  referred  to  today  as  the  asymptotic,  threshold/ambiguity, 
and  no  information  regions. 


can  interpret  the  range  estimate  as  a  one-dimensional  output  of  the  demodulation/decoding  of  a 
higher  dimensional  input  (the  waveform). 

The  interpretation  of  the  threshold  effect  in  DOA  estimation  error  from  the  perspective  of 
Shannon’s  thesis  is  rather  transparent  in  light  of  Woodward’s  work,  but  noteworthy  nonetheless. 
The  relative  phase  from  sensor- to-sensor  of  the  signal  component  (known  as  the  signal  steering  vec¬ 
tor,  array  response  vector,  or  system  function  vector,  etc.)  can  be  interpreted  as  a  multidimensional 
encoding  of  a  one-dimensional  message  (the  signal’s  arrival  angle).  The  same  interpretation  applies 
to  matched  field  processing,  for  example,  where  the  low-dimensional  encoded  message  would  be  the 
source  range  and  depth. 

Van  Trees  [59]  discusses  the  threshold  effect  in  the  context  of  classical  estimation  theory. 
This  discussion  is  embedded  primarily  in  an  example  problem  (Example  2  of  Section  4.2.3  Non¬ 
linear  Estimation),  but  provides  a  general  strategy  for  predicting  MSE  performance  of  nonlinear 
estimation  schemes  when  local  bounds  like  the  CRB  are  no  longer  useful.  The  example  problem 
deals  with  ML  estimation  of  waveform  parameters  where  the  waveforms  are  those  of  a  nonlinear 
pulse  frequency  modulation  scheme,  and  the  general  strategy  proposed  for  threshold  region  MSE 
prediction  is  referred  to  herein  as  MIE.  A  detailed  description  of  MIE  shall  be  deferred  to  Section 
4  of  this  report,  but  the  short  description  is  that  it  builds  on  a  judicious  decomposition  of  the 
integral  for  the  MSE  of  the  parameter  estimate.  Van  Trees  approximated  the  necessary  interval 
error  probabilities  with  standard  asymptotics  and  use  of  well  known  bounds  on  the  error  function. 

Rife  and  Boorstyn  (RB)  [43]  build  on  Van  Trees’  approach  to  predict  the  threshold  region 
performance  of  a  Discrete  Fourier  Transform  (DFT)-based  approximate  implementation  of  an  ML 
estimator.  RB  considers  the  estimation  of  the  frequency  of  a  single  pure  sinusoid  in  additive  white 
Gaussian  noise  (AWGN)  where  the  sinusoid  is  assumed  deterministic  (and  therefore  represented 
in  the  mean  of  the  distribution).  True  ML  requires  a  full  search  over  a  continuous  interval  of 
frequencies.  RB  considers  a  DFT-based  approach  that  restricts  the  ML  continuous  search  to  a 
finite  set  of  discrete  sample  points.  When  sampled  at  Nvquist,  the  DFT  samples  the  mainlobe  of 
the  ambiguity  function  (a  sine  in  this  case),  and  the  zeros  of  the  ambiguity  function  in  the  sidelobe 
region.  Thus,  RB  considers  the  impact  of  restricting  the  ML  search  to  discrete  sampling  of  the 
mainlobe  and  the  zeros  of  the  ambiguity  function.  Since  the  DFT  bins  corresponding  to  the  zeros 
of  the  ambiguity  function  are  orthogonal  to  the  mainlobe,  the  analysis  of  RB  (derivation  of  the 
interval  error  probabilities,  or  outlier  probabilities  as  RB  refers  to  them)  simplifies  significantly. 
This  landmark  paper  popularized  Van  Trees’  MIE  approach  to  below  threshold  MSE  prediction. 

Steinhardt  and  Bretherton  [48]  took  the  analysis  of  RB  two  steps  further  by  (i)  simplifying 
the  expressions  for  the  error  probabilities,  and  (ii)  providing  for  what  appears  to  be  the  first  time 
a  method  for  computing  directly  in  closed  form  the  threshold  SNR  of  the  RB  analysis.  Direct 
calculation  of  the  threshold  SNR  simply  means  that  an  initial  prediction  of  the  MSE  curve  as  a 
function  of  SNR  (from  which  the  threshold  SNR  is  easily  determined)  is  bypassed. 

Tufts,  Kot,  and  Vaccaro  (TKV)  [54,  55]  developed  a  procedure  for  computing  the  threshold 
SNRs  for  the  estimation  of  the  frequencies  of  multiple  sinusoids  by  linear  prediction.  Their  proce¬ 
dure  can  be  viewed  as  a  variant  of  MIE.  This  creative  analysis  appears  to  be  the  first  to  recognize 
the  general  applicability  of  MIE  to  a  large  class  of  nonlinear  estimation  schemes.  It  is  also  the  first 


6 


to  jointly  quantify  the  below  threshold  SNR  performance  of  the  estimation  of  multiple  sinusoids.  As 
mentioned  in  the  introduction,  application  of  MIE  requires  good  approximation  of  two  quantities: 
(i)  the  asymptotic  MSE  of  the  algorithm,  and  (ii)  the  interval  error  probabilities.  TKV  obtained 
the  asymptotic  MSE  of  linear  prediction  by  applying  a  matrix  perturbation  analysis  to  essentially 
approximate  the  Taylor’s  series  of  the  OSF.  The  linear  and  quadratic  terms  are  retained  to  com¬ 
pute  the  desired  local  error  MSE  of  the  parameter  estimates.  TKV  approximated  the  required  error 
probabilities  by  first  defining  a  useful  upperbound  on  the  desired  probability  that  is  much  easier  to 
work  with  analytically,  and  secondly  assuming  the  essential  random  variables  involved  in  the  error 
event  are  well  approximated  as  Gaussian.  The  latter  assumption  allows  use  of  the  well  tabulated 
error  function.  The  procedure  defined  by  TKV  is  classified  herein  as  a  variant  of  MIE  because  no 
semblance  of  intervals  of  error  are  ever  defined.  The  recognition  of  the  decomposition  of  the  total 
MSE  as  the  sum  of  two  contributions  (global  errors  and  local  errors),  as  given  in  equation  (127)  on 
p.  282  of  [59],  but  repeated  in  a  form  similar  to  that  used  in  RB  [43],  however,  is  present  in  their 
work. 


An  interesting  analysis  related  to  TKV  is  the  work  of  Kaveh  and  Wang  (KW)  [25]  that  focuses 
on  the  low  SNR  performance  of  the  Multiple  Signal  Classification  (MUSIC)  algorithm  [4,  52,  61] 
and  the  minimum-norm  (min-norm)  algorithm.  These  are  subspace  based  algorithms  that  exploit 
the  singular  value  decomposition  (SVD)  as  does  the  linear  prediction  algorithm  examined  by  TKV. 
KW  focuses  on  quantifying  how  well  both  algorithms  detect  (estimate  the  number  of  planewaves), 
estimate  the  signal  subspace,  and  resolve  closely  spaced  sources.  KW  approximates  the  probability 
of  resolution  and  uses  this  as  a  lower  bound  on  the  actual  threshold  SNR  for  closely  spaced  sources. 
KW  primarily  builds  upon  well  known  asymptotic  results  on  the  statistics  of  the  eigenvectors  and 
eigenvalues  of  a  complex  Wishart  random  matrix  to  quantify  their  desired  performance  measures. 

The  recent  work  of  Xu  et  al.  [67]  [71]  that  considers  source  localization  error  performance  of 
MFP  methods  in  underwater  acoustic  environments,  the  work  of  Athley  [2,  3]  that  considers  the 
ML  DOA  estimation  of  multiple  planewave  signals  in  white  noise,  and  that  of  Boyer  et  al.  [12,  13] 
that  also  considers  the  ML  DOA  estimation  of  planewaves  in  white  noise,  extend  the  work  of  [43]  in 
several  ways.  First,  the  restriction  that  the  ML  search  be  limited  to  DFT  bins  sampled  at  Nyquist 
is  relaxed.  This  leads  to  performance  that  better  approximates  true  ML,  but  it  likewise  complicates 
the  calculation  of  the  interval  error  probabilities  because  the  sidelobes  of  the  ambiguity  function 
and  the  mainlobe  are  statistically  dependent  in  general;  second,  the  signal  model  is  no  longer 
restricted  to  be  a  sinusoid  (or  planewave);  and  third,  analysis  was  extended  to  include  stochastic 
signal  models  (signal  component  represented  in  the  covariance  of  the  distribution).  Richmond  [37] 
has  recently  extended  these  analyses  to  ML  signal  parameter  estimation  within  an  adaptive  array 
context. 

Before  closing  this  section,  it  is  worthwhile  to  point  out  that  researchers  have  considered 
Bayesian  bounds  to  address  the  contribution  of  global  errors  to  the  overall  MSE  of  ML.  Bounds 
receiving  the  most  attention  include  the  Ziv-Zakai  bound  (ZZB),  Bayesian  CRB  [59],  Barankin 
bound  [36],  and  the  Weiss- Weinstein  bound.  Recent  work  has  extended  the  theory  of  ZZBs  to 
include  vector  parameters  [8,  9].  The  Bayesian  framework  is  attractive  as  its  applicability  is  not 
limited  to  unbiased  estimators,  and  it  provides  a  viable  means  of  including  contributions  of  global 
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errors  to  the  overall  MSE3.  Loosely  speaking,  the  nature  of  the  ML  search  is  effectively  mimicked 
by  these  bounds  by  randomizing  the  signal  parameter  of  interest,  e.g the  signal’s  angle  of  arrival 
becomes  a  random  variable  with  a  known  distribution.  This  is  done,  however,  at  the  expense 
of  the  classical  assumption  of  ML  theory  that  the  unknown  parameters  are  deterministic.  The 
ability  to  provide  system  performance  prediction  for  a  specific  parameter  value  or  given  set  of 
parameter  values  is  a  desired  characteristic  of  a  performance  bound.  Indeed,  the  non- Bayesian 
CRB  has  enjoyed  much  popularity  in  practice  as  a  system  design  tool,  in  part  because  it  provides 
such  flexibility.  This  flexibility  is  not  offered  by  Bayesian  bounds  by  construction.  In  addition, 
Bayesian  bounds  are  often  computationally  expensive,  do  not  easily  support  inclusion  of  the  effects 
of  estimating  nuisance  parameters  (like  an  unknown  data  covariance  matrix),  and  at  best  predict 
the  performance  of  ML  estimation.  MIE  in  contrast  is  an  approach  to  MSE  prediction  that  is 
algorithm  specific,  maintains  the  assumption  that  parameters  are  deterministic  yet  unknown,  and 
can  support  the  presence  of  nuisance  parameter,  although  the  analyses  can  be  nontrivial. 


^Extensions  of  classical  local  bounds  to  biased  estimators  do  exist.  Their  use,  however,  is  strictly  limited  in  that 
good  knowledge  of  the  bias  is  required  in  order  to  compute  the  bound  [59]. 
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3.  THE  CAPON  AND  BARTLETT  ALGORITHMS  FOR  DOA 

ESTIMATION 


3.1  The  Spectral  Estimation  Problem 

The  Capon  and  Bartlett  algorithms  are  well  known,  and  their  performance  has  been  studied 
extensively.  The  goal  of  this  section  is  to  establish  notation  and  to  briefly  summarize  the  algorithms 
and  known  results  on  their  large  sample  local  error  MSE  performance. 

The  fundamental  problem  of  spectral  estimation  can  be  stated  as  follows:  Given  a  finite 
length  data  record,  ie.,  a  sample  function  of  a  random  process  that  is  at  least  wide  sense  station¬ 
ary  observed  over  some  finite  interval  of  time,  estimate  the  power  distribution  over  frequency.  The 
Capon  and  Bartlett  algorithms  represent  filterbank  approaches  to  this  classic  problem.  Although 
originally  posed  for  time  series  analysis  where  frequency  refers  to  temporal  frequency,  the  power 
spectrum  of  interest  for  the  DOA  estimation  problem  is  with  respect  to  spatial  frequency  or  angle. 
To  this  end,  data  is  obtained  from  an  array  of  N  sensors  distributed  in  space.  Each  data  observa¬ 
tion  taken  across  the  array  (called  a  spatial  snapshot)  is  modeled  as  an  N  x  1  zero  mean  complex 
Gaussian  vector4  with  representation  x  =  Sv(6)  +  n  where  the  array  response  of  the  signal  of 
interest  associated  with  DOA  parameter  6  is  given  by  v(0),  and  its  complex  amplitude  is  Gaussian 
distributed  such  that  S  ~  CA/’i(0,  a|),  and  the  colored  noise5 *  is  denoted  by  n  having  covariance 
Ejnn^}  =  R.y,  where  E{-}  denotes  the  statistical  expectation.  Note  that  £{x}  =  0  and  that 
E{xxh}  =  R  =  Ryv  T  (t^v(6)vh  (6).  A  finite  set  of  L  array  observations  is  accrued  over  time  and 
assembled  in  a  data  matrix:  X  =  [x(l)|x(2)|  •  •  •  |x(L)],  where  x(Z)  ~  CM  n{  O.R),  /  =  1,2,  ...,L. 
These  snapshots  often  represent  measurements  of  the  narrowband  spatial  complex  envelope  whose 
real  and  imaginary  parts  are  respectively  composed  of  the  in-phase  and  quadrature  component  of 
demodulated  data.  The  spatial  snapshots  are  used  to  form  the  unnormalized  data  spatial  covari¬ 
ance  estimate  R  =  XX 11  from  which  the  Capon  and  Bartlett  algorithms  generate  power  spectral 
estimates. 

The  Capon  and  Bartlett  DOA  angle  estimates  are  obtained  as  the  arguments  of  the  largest 
peaks  of  the  estimated  spatial  power  spectra.  If  P(0)  represents  the  estimated  spectrum  as  a 
function  of  angle,  then  the  maximum  output  provides  an  estimate  of  the  signal  power  cr|,  and  the 
signal  DOA  estimate  is  given  by  the  scan  value  of  0  that  achieves  this  maximum;  namely, 

6  =  arg  max  P{9)  ( 1 ) 

6 

(assuming  a  single  signal  is  present).  It  shall  be  assumed  that  K  signals  are  present  in  the  data,  and 
that  the  Capon/Bartlett  parameter  estimates  0*,  k  =  1,2,...,  K  are  obtained  as  the  arguments  of 

4The  notational  convention  adopted  is  as  follows:  italics  indicates  a  scalar  quantity,  as  in  A\  lower  case  boldface 
indicates  a  vector  quantity,  as  in  a;  upper  case  boldface  indicates  a  matrix  quantity,  as  in  A.  The  n-th  row  and  m- th 
column  of  matrix  A  will  be  indicated  bv  [A]n,m.  Variables  will  be  assumed  complex  in  general,  but  some  will  be  real 
(obvious  from  context).  Re(A)  is  the  real  part  of  A  and  Im(A)  is  the  imaginary  part.  The  complex  conjugation  of 
a  quantity  is  indicated  by  a  superscript  *  as  in  A*.  The  matrix  transpose  is  indicated  by  a  superscript  T  as  in  A1  , 
and  the  complex  conjugate  plus  matrix  transpose  is  indicated  by  a  superscript  H  as  in  Al!  =  (A7  )*. 

5It  is  assumed  that  the  colored  noise  does  not  contain  any  low  rank  directional  signals  in  the  search  domain,  since 

the  total  number  of  signals  in  search  space  is  assumed  known. 
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the  K  largest  peaks  of  P{0). 


3.2  The  Bartlett  Algorithm:  A  Conventional  Beamforming  Approach 

Bartlett  recognized  the  inherent  tradeoff  between  resolution  and  variability  in  the  periodogram 
(Fourier  Transform  of  a  sample  data  record)  spectral  estimate  and  recommended  averaging  over 
shorter  time  Fourier  Transforms  as  a  means  of  reducing  the  high  variance  at  the  expense  of  some 
resolving  capacity.  The  Bartlett  algorithm  estimates  power  spectra  over  angle  by  averaging  the 
power  out  of  a  bank  of  conventional  spatial  filters  (beamformers).  Each  conventional  spatial  filter 
is  a  crude  narrowband  (in  spatial  frequency)  filter  tuned  to  and  centered  on  a  specific  angle  of 
interest.  The  output  of  each  filter  is  power  averaged  over  time,  and  with  proper  normalization 
provides  an  estimate  of  the  true  underlying  power  spectral  density  such  that  when  integrated  over 
the  spatial  frequency  band  of  interest  its  area  provides  an  estimate  of  the  power  of  the  random 
process  contained  in  that  band.  The  resulting  spectra  plotted  as  a  function  of  the  filter  center 
frequency  produces  what  is  classically  referred  to  as  a  smoothed  periodogram. 

Specifically,  let  the  conventional  beamforming  weight  steered  to  an^le  0  for  an  array  with 
element  spatial  locations  zn,  n  =  1, 2, . . . ,  N  be  given  by  v(Q)  =  [ejk^Zl,  22 , . . . ,  ejk" Zv] 7 , where 
ktf  =  (27r/A)a(0)  is  the  wavenumber  vector.  a(0)  is  the  3x1  unit  vector  pointing  in  the  assumed 
direction  of  field  propagation,  A  =  c/f  is  the  wavelength  at  temporal  frequency  /,  and  c  is  the 
wave  propagation  speed.  Note  that  for  a  uniform  linear  array  (ULA),  this  spatial  filter  weight 
can  be  easily  implemented  with  the  Fast  Fourier  Transform  (FFT).  The  Bartlett  spectral  estimate 
evaluated  at  spatial  frequency  (or  angle)  0  is  given  by 

1  L  1 

PBartlett(O)  =  T  E  lV"(0)X(O|  =  T  '  v"(0)R v(0)  (2) 

L  /  =  1 

and  its  ambiguity  function  is  defined  as  ipsa 


3.3  The  Capon  Algorithm:  An  Adaptive  Beamforming  Approach 

Capon  recognized  that  conventional  filterbank  techniques  for  spectral  estimation  rely  too 
heavily  upon  the  passband  and  stopband  properties  of  a  Fourier  Transform  frequency  bin  filter: 
namely,  the  FFT  bin  has  a  frequency  response  given  by  a  sine  function.  Taper  design  has  been 
exploited  to  reduce  sidelobes  (and  leakage),  at  the  expense  of  a  wider  mainlobe  (and  reduced 
resolution)  [29,  52].  When  multiple  signals  are  present,  sidelobe  leakage  and  mutual  signal  interfer¬ 
ence  can  introduce  biases  in  DOA  estimates  obtained  via  conventional  techniques  like  the  Bartlett 
algorithm  (see  Section  9).  Capon  proposed  designing  these  filters  optimally;  namely,  a  linear  nar¬ 
rowband  filter  designed  to  pass  the  desired  signal  undistorted,  while  minimizing  the  power  from 
all  other  spatial  frequency  bands  (and  thus,  other  sources  of  interference).  Formally,  the  following 
constrained  optimization  problem  for  the  filter  weight  w  was  suggested 

R,- 1  0) 

min  w^Rw  such  that  w"v(0)  =  1  =>  w mvdr  =  i  (3) 

w  vrt  (t/)R-1  v(t/) 
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that  yields  the  well  known  MVDR  filter  as  its  solution.  Note  that  the  optimal  filter  weight  depends 
on  the  data  covariance  matrix  R.  This  filter  is  often  interpreted  in  the  classical  sense  as  a  data- 
dependent  taper  [15,  17].  By  construction,  this  optimal  filter  will  cancel  all  spatially  coherent 
energy  from  all  directions  other  than  the  scan  direction  0.  The  filter  output  w  my  dr*  likewise 
provides  the  ML  estimate  of  the  complex  signal  amplitude  S  [33,  49,  39,  31].  The  average  output 
power  is  given  by 


(I 


E \  |  WMVDRX 


vH(0)R-1v(0) 


—  H’Capon 


(9)  = 


1 


VW(0)RyV(0) 


+  °s 


(4) 


where  the  Capon  ambiguity  function  ipcapon(Q)  has  been  defined.  The  last  equality  in  (4)  holds 
only  when  the  signal  array  response  in  R  perfectly  matches  that  used  to  form  the  weight  vector 
wmvdr,  and  when  R  is  perfectly  known.  Capon,  therefore,  reasoned  that  for  large  enough  sample 
support,  an  estimate  of  the  covariance  R  can  be  used  with  (4)  to  estimate  the  signal  power  <r| 
and  the  corresponding  signal  parameter  0  (or  set  of  signal  parameters  6  =  [#ij,  #1,2, ....  0\,aY  as 
in  [15,  21,  24]).  Using  the  covariance  estimate  R  =  XX/;,  Capon  proposed  the  following  power 
spectral  estimator 

Pc°’™m  =  L-N  +  1  '  <5) 

Pcapan{0)  can  be  further  normalized  to  ensure  that  it  is  a  true  power  spectral  density  [29,  52]. 


3.4  Asymptotic  Local  Error  MSE  Performance 

Consider  the  case  of  a  single  planewave  signal  in  AWGN  and  note  that  at  high  enough 
SNRs,  namely,  those  above  the  threshold  SNR,  the  angle  search  indicated  in  (1)  will  result  in  a 
DOA  estimate  obtained  from  the  local  neighborhood  of  the  true  global  maximum  of  the  ambiguity 
function  with  near  probability  1.  At  such  SNRs  it  is  possible  to  characterize  the  average  error  or 
jitter  in  DOA  estimates  about  the  true  value  by  exploiting  Taylor’s  theorem.  For  ML  estimation 
of  DOAs,  this  jitter  error  is  characterized  by  the  CRB  (it  will  be  demonstrated  in  Section  9  that 
the  Bartlett  algorithm  in  fact  produces  the  ML  DOA  estimate  for  the  single  signal  case  in  white 
noise).  The  same  ideas  hold  in  the  multiple  signal  case,  where  the  estimate  of  the  k- th  signal  angle 
will  be  obtained  from  the  local  neighborhood  of  the  fc-th  largest  peak  of  the  ambiguity  function 
with  near  probability  1  (assuming  all  signals  have  large  SNRs).  The  following  sections  describe  the 
fundamental  results  of  analyses  that  compute  this  jitter  error,  ie.,  the  asymptotic  local  error  MSE 
of  the  DOA  estimate  of  the  k- th  signal. 

It  will  be  assumed  in  this  section  that  for  well  separated  (by  at  least  a  beamwidth)  signals, 
SNRs  exceed  estimation  threshold,  and  that  for  closely  spaced  signals,  SNRs  are  large  enough  for 
all  signals  to  be  mutually  resolved  with  zero  probability  of  intersource  errors.  In  addition,  it  is 
assumed  that  no  pairs  of  perfectly  coherent  signals  are  present  (coherent  sources  are  not  resolvable 
with  the  Capon  algorithm  [61]),  and  that  the  total  number  of  signals  present  in  the  data  is  known0. 

f>Subspace  rank  (model  order)  determination  is  an  important  topic  and  active  area  of  research,  but  will  not  be 
pursued  herein. 
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3.4.1  Bartlett  Algorithm 

Hawkes  and  Nehorai  (HN)  [21]  developed  the  theory  based  on  Taylor’s  theorem  for  predicting 
the  asymptotic  local  error  MSE  performance  of  DOA  estimates  derived  from  the  Bartlett  spectral 
estimator.  To  summarize  the  results  of  HN,  let  6^,  k  =  1,2,...,  K  represent  the  true  DOA  angle 

values  of  the  K  signals  present  in  the  data.  Let  k  =  1,2 . K  represent  the  asymptotic 

(L  — ►  oo)  estimates  of  these  K  angles,  i.e.,  the  locations  of  the  K  largest  peaks  of  the  ambiguity 
function  ipBartlett{0 )  (that  are  not  necessarily  equal  to  6^  for  all  SNR).  The  MSE  is  the  sum  of  the 
bias  squared  and  variance.  Thus,  the  MSE  of  the  estimate  0 ^  is  given  by 

£{(&-0fc)2}  =  [^{(^.-^)}]2  +  £:|(^-4)2}  (6) 

where  the  bias  can  be  written 

E{dk-ek)  =  E  {4  -  §k }  +  (ek  -  ek) .  (7) 

HN  derived  expressions  for  the  asymptotic  bias  and  variance  needed  to  approximate  the  MSE  in 
(6).  Denote  the  first  three  derivatives  of  the  array  response  with  respect  to  the  scan  variable  0, 
respectively,  as 


d'2v(8) 


v  (9)  = 


d'W(9) 

d(P 


(8) 


As  mentioned  in  [21],  the  additional  bias  term  for  the  Bartlett  DOA  estimate  is  negligible  in  most 
cases,  i.e.,  E{8k  —  9k}  ~  0.  When  the  asymptotic  bias  is  small  and  nonzero,  it  can  be  approximated 
via 


(4  -  Ok) 


-Re  [v"(gfc)Rv(flfc)] 


Re  [v"(0fc)Rv(0fc)  +  v"(0*)Rv(0*)] ' 


(9) 


If  the  bias  is  large,  then  (9)  can  be  grossly  inaccurate.  The  bias,  however,  can  always  be  easily 
inferred  directly  from  the  location  of  the  true  k-th  peak  of  the  ambiguity  function  Bartlett {0), 
relative  to  the  true  A>th  signal  angle.  The  asymptotic  variance  is  given  by 


E{(*-*)’}  “51- 


vH (9k)Rv(0k)vH (9k)Rv(9k)  -  i," (9k)Rv(9k)vH (0k)Rv(9k) 


II  i 


.Hi 


Re 


v"(9k)Rv(9k)  +  vH  (9k)Rv(9k)  1 


(10) 


The  total  asymptotic  local  error  MSE  approximation  obtained  for  the  Bartlett  DOA  estimates 
using  (6)-(10)  shall  be  herein  denoted  by  cr2HN(0^). 


3.4.2  Capon  Algorithm 

The  large  sample  (L  N)  asymptotic  local  error  MSE  performance  of  the  Capon  signal 
parameter  estimator  has  been  theoretically  analyzed  by  several  authors.  Stoica  et  al.  [51],  VB  [57], 
and  HN  [21]  exploit  Taylor’s  theorem  and  complex  gradient  methods  to  approximate  the  MSE. 
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VB  provides  an  additional  bias  term  via  a  second  order  Taylor  series  expansion  that  is  intended  to 
further  capture  finite  sample  effects  and  therefore  a  broader  range  of  values  for  L  and  SNR. 

To  summarize  the  results  of  VB,  let  0 k  =  1,2,...,  K  represent  the  true  parameter  values 
of  the  K  signals  present  in  the  data.  Let  6 k  =  1,2,...,  IV  represent  the  asymptotic  ( L  — ►  oo) 
estimates  of  these  K  parameters,  i.e .,  the  locations  of  the  K  largest  peaks  of  the  ambiguity  function 
Caponi &)  (that  are  not  necessarily  equal  to  6 p  for  all  SNR).  The  MSE  of  the  estimate  Op  is  given 
in  form  by  (6)  with  bias  given  by  (7).  Recalling  the  derivatives  given  in  (8),  define  the  function 
fcapon{0 )  and  denote  its  second  and  third  order  derivatives,  respectively,  as 

f Capon  =  VW(0)R->V(0)  =  1  /4>Capon(0) 

fcapon(0)  =  2Re  [v"(6)R-lx(9)  +  v"(0)R -*v((?)]  (11) 

f Capon  (0)  =  2Re  [3v*(0)R ~'v(9)  +  vw(«)R->  V  (0)] 

and  define  the  matrix  B (0*.)  as 

b(4)  =  v(4)v"(4)  +  v(4)vw(4)-  (12) 


The  following  expressions  for  the  bias  and  asymptotic  variance  of  the  Capon  algorithm  signal 
parameter  estimates  were  derived  by  VB  [57]: 


-Rc 


E  {  (4  -  4)  }  * 


"(^)R-!v(^)  +  v"(0fc)R-1v(^.)] 

vw(4)R-,B(4)R"1v(4) 


2{L  —  N) 


T  Re 


{L-N -\)(L-N +  \)m 

f Capon  i^k) 


[fCapon(Ok)}2 


2  f Capon  i@k) 


E{(ok-oky} 


(13) 


ft*  0VM  2{L~N)  ] 

Re 

vh(4)r-1b(4)r-1v(4) 

\Vfc  9k)  J  _  [(L-iV-l)(L-yV  +  l)_ 

[fcapon{9k)Y 

(14) 


The  Capon  large  sample  MSE  is  obtained  by  using  relations  (11)  (14)  in  the  expression  for  MSE 
(6)  (7).  This  MSE  approximation7  shall  be  herein  denoted  by  the  symbol  o-‘yB(0p). 

Note  in  general  that  both  MSE  approximations,  &2HN{0k)  for  the  Bartlett  algorithm  and 
cryB(0k)  for  the  Capon  algorithm,  are  functions  of  the  array  geometry,  number  of  sensors  V, 
sample  support  L,  true  data  covariance  R  (includes  colored  noise  and  true  array  manifold  of 
signals),  assumed  scanning  vectors  v(0),  and  SNRs.  Thus,  the  asymptotic  MSE  performance  can 
be  explored  versus  combination  of  these  parameters. 

It  was  observed  in  this  analysis  that  the  additional  bias  term  was  often  negligible  (as  mentioned  in  [21]).  It  was 
likewise  found  that  the  approximate  bias  expression  given  by  the  first  equation  in  (13)  can  be  grossly  inaccurate  when 
the  bias  is  large.  A  reliable  estimate  of  the  bias  (Ok  —  Ok),  however,  is  easily  inferred  directly  from  the  location  of  the 
peak  of  the  ambiguity  function  ^ caPon(6 )  relative  to  the  true  signal  parameter  value  Ok- 
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4.  AN  INTERVAL  ERROR-BASED  METHOD  OF  THRESHOLD  REGION 
MEAN  SQUARED  ERROR  PREDICTION 


This  section  describes  the  method  of  interval  errors  (MIE)  for  MSE  prediction  and  its  adapta¬ 
tion  to  the  Capon  (and  Bartlett)  algorithm.  Application  of  MIE  to  the  single  signal  case  (K  =  1)  is 
restricted  to  estimation  scenarios  that  are  at  least  asymptotically  (SNR— ►  oo)  well-posed  estimation 
problems,  consisting  of  ambiguity  functions  that  possess  a  unique  global  maximum  over  the  signal 
parameter  search  space.  While  this  technique  can  be  used  for  arbitrary  array  configurations  and 
signal  models  (more  general  than  planewaves),  it  cannot  predict  performance  when  the  ambiguity 
function  does  not  possess  a  unique  global  maximum  over  the  parameter  search  domain  at  least  for 
very  large  SNRs.  Extension  of  MIE  to  the  multisignal  case  (A'  >  1)  requires  that  the  ambiguity 
function  possess  K  unique  local  maxima  due  to  the  K  sources. 

The  following  discussion  assumes  that  no  signal  modeling  errors  are  present.  Mismatch  will  be 
considered  in  Section  8.  MIE  is  first  described  for  the  case  of  a  single  source  present.  Its  extension 
to  multiple  sources  is  described  subsequently.  The  steps  for  computing  the  required  pairwise  error 
probabilities  for  both  the  Capon  estimator  and  the  Bartlett  are  then  summarized.  Simple  large 
sample  approximations  of  these  pairwise  error  probabilities  based  on  the  complementary  error 
function  are  derived  in  Section  4.2.1. 

The  following  discussion  focuses  on  the  Capon  algorithm.  The  parallel  application  to  the 
Bartlett  algorithm  is  transparent;  however,  a  summary  of  the  analogous  results  is  provided  at  this 
section’s  end. 


4.1  Method  of  Interval  Errors 

MIE  builds  upon  the  two  regions  of  the  composite  MSE  curve  of  Figure  1  that  are  given  by 
the  asymptotes  of  the  SNR;  namely,  the  no  information  (SNR— »  0)  and  asymptotic  (SNR— ►  oo) 
regions.  Define  the  conditioning  event 

A  =  {True  source  parameters  are  0&,  k  =  1, 2 . A"}  .  (15) 

MIE  decomposes  the  MSE  expression  into  two  components:  “no  interval  errors”  (NIE)  and  “interval 
errors”  (IE) 


E  ( (ek  -  ^)1-4  =  [  ps(  h  =  do  .4)  (do  -  dkf  ddo 

r vi  \  ( vi  i  (16) 

=  Pr  ( NIE  1 4)  £  <  \0k-9k  J  NIE.  .4  >  +  Pr  ( IE  1 4)  £  <  \9k-9k)  IE..O 


(see  equation  (127)  on  p.  282  of  [59]).  The  parameter  search  space,  z.e.,  the  scanning  domain  for 
0,  is  divided  into  disjoint  mutually  exclusive  intervals  based  on  the  characteristics  of  underlying 
ambiguity  function  ^Capon{Q)  that  innately  drives  the  character  of  the  stochastic  process  generated 
by  the  OSF  Pcapon{0 )  011  each  scan.  Note  from  (4)  that  'ipCapon{@ )  depends  on  R,  and  is  a  function 
of  the  K  SNRs  of  the  K  signals  present.  Consequently,  this  SNR  dependence  can  affect  the  interval 
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choice  and  must  be  accounted  for  in  general.  When  the  SNR  dependence  is  weak  however,  then 
a  single  choice  of  intervals  will  suffice  for  all  SNRs  with  no  need  to  adjust.  For  the  presentation 
that  follows,  weak  dependence  is  assumed.8  The  ambiguity  function  used  for  interval  selection  is 
obtained  from  an  arbitrary  choice  of  some  high  value  of  SNR. 

4.1.1  Single  Source,  A' ’  =  1 

A  large  enough  SNR  (how  large  is  scenario  dependent)  leads  to  an  ambiguity  function  with  a 
dominant  peak  located  at  the  true  parameter  value,  say  0\.  and  several  false  peaks  located  elsewhere 
within  the  domain  of  interest.  An  example  ambiguity  function  is  illustrated  in  Figure  2.  MIE  divides 
the  parameter  search  space  into  disjoint  neighborhoods  or  intervals  around  all  local  maxima  of  the 
ambiguity  function  r/j Capon{0 ).  and  interprets  the  resulting  MSE  in  terms  of  these  intervals.  When 
signal  parameter  estimates  fall  inside  the  interval  containing  the  global  maximum  0\,  then  it  is 
said  that  “no  interval  errors”  have  occurred.  Such  estimates  result  in  the  smallest  contributions 
to  the  MSE,  often  referred  to  as  local  errors  (jitter  errors).  This  NIE  component  represents  the 
contribution  dominating  the  MSE  within  the  asymptotically  high  SNR  regime.  When  considering 
ML  estimation,  for  example,  the  CRB  is  often  a  good  predictor  of  the  ML  MSE  performance  in  this 
region,  describing  the  small  perturbations  or  jitter  in  the  resulting  signal  parameter  estimate  as  the 
ML  search  resides  with  near  probability  1  in  a  small  neighborhood  of  the  true  peak  of  the  ambiguity 
function.  The  large  sample  MSE  approximation  obtained  via  (JyB(0\)  (as  defined  in  Section  3.f.2) 
for  signal  parameter  estimates  derived  from  the  Capon  algorithm  will  be  used  to  describe  this  NIE 
component  contribution  to  the  overall  MSE.  This  is  appropriate  by  virtue  of  Taylor’s  theorem  upon 
which  cryB(6\)  is  based.  The  truncated  Taylor  series  approximates  the  local  behavior  of  the  Capon 
OSF  about  the  point  of  expansion  6\. 

When  signal  parameter  estimates  fall  outside  the  interval  consisting  of  the  local  neighborhood 
around  the  global  maximum  of  4> Caponi Q)  located  at  0 i,  then  it  is  said  that  an  “interval  error”  has 
occurred.  Such  estimates  result  in  the  largest  contributions  to  the  MSE,  contributing  what  are 
referred  to  as  global  errors.  This  IE  component  results  as  a  consequence  of  the  Capon  algorithm 
choosing  maxima  of  the  OSF  that  can  be  traced  to  false  peaks  (those  local  maxima  not  corre¬ 
sponding  to  0\)  of  the  underlying  ambiguity  function.  This  component  dominates  in  the  low  to 
very  low  SNR  regimes  and  must  be  approximated  numerically.  An  adequate  approximation  can  be 
obtained  by  exploiting  a  property  that  is  characteristic  of  nonlinear  estimation  schemes  consisting 
of  an  OSF  driven  by  an  underlying  multimodal  ambiguity  function;  namely,  that  the  pdf  of  the 
resulting  signal  parameter  estimate  tends  to  aggregate  its  density  in  parameter  space  0  around 
the  local  maxima  of  the  multimodal  ambiguity  function,  especially  at  SNRs  in  the  vicinity  of  the 
estimation  threshold  region  (see  histogram  in  Figure  2  that  corresponds  to  the  ambiguity  function 
shown  to  its  left;  also,  see  analogous  histograms  for  MUSIC  in  Figure  9.30,  p.  1199  of  [hi],  for 
example).  Consequently,  the  continuous  pdf  of  the  Capon  signal  parameter  estimate  can  be  well 
approximated  by  a  discrete  probability  mass  function  (pmf)  (or  a  finite  set  of  Dirac  delta  functions 
in  continuous  parameter  space),  where  the  masses  (or  delta  function  areas)  are  given  by  interval 
error  probabilities.  Specifically,  let  all  local  maxima  of  the  ambiguity  function  4’CaponiQ)  within  the 


HThis  assumption  is  well  satisfied  for  well  separated  signals  in  white  noise.  Colored  noise  and/or  closely  spaced 
signals  can  require  that  this  dependence  be  accounted  for  (see  [30],  for  example). 
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CAPON  AMBIGUITY  FUNCTION  HISTOGRAM  OF  CAPON  ALGORITHM 
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Figure  2.  Example  Capon  ambiguity  function  and  histogram  with  true  signal  at  0 1  =  90  degrees. 


signal  parameter  domain  of  interest  when  evaluated  at  a  sufficiently  large  SNR  (large  enough  such 
that  the  true  peak  is  evident  and  located  at  6\)  be  given  by  the  finite  set  M  \  =  {  0  \  Q\,  62, . . . ,  6  m  }, 
where  0\  is  the  global  maximum  location  and  6k  for  k  =  2,3,  ...,71/  are  all  other  local  maxima; 
with  this  discretization  one  can  effectively  approximate  the  continuous  scanning  search  over  an  in¬ 
finite  number  of  possible  outcomes  with  a  A/-ary  hypothesis  testing  over  a  finite  number  of  possible 
estimates.  The  total  MSE  for  this  Capon  parameter  estimate  can  be  approximated  by9 * 11 


£<[(0.  _6,|)'|'4}  -  [1_  X>(^1  =  ^"'l'4) 

A/ 

+  £p(*l  =0m|.4)  (Om-Oi)2 


<Tvb(9  1) 


777—2 


(17) 


where  the  NIE  component  is  described  by  the  results  of  VB,  and  the  IE  component  has  been 
discretized.  The  interval  error  probability  p  (^0 \  =  6m  represents  the  likelihood  of  the  Capon 


search  algorithm  choosing  the  false  peak  located  at  6  =  6m  as  an  estimate,  when  the  true  signal  is 
located  at  parameter  value  0  =  9 Technically  speaking,  the  interval  error  probability  is  obtained 
by  integrating  the  pdf  of  the  parameter  estimate  over  the  error  interval.  Since  the  pdf  is  unknown, 
the  approximation  in  (17)  focuses  on  the  single  point  in  each  interval  most  likely  to  contend  with 
the  global  maximum  at  6\\  namely,  the  local  maximum  contained  in  each  interval.  This  single  point 


9 Although  negligible  in  the  threshold  region,  the  following  adjustment  to  the  local  error  contribu¬ 

tion  is  often  necessary  in  the  no  information  region  due  to  the  approximate  nature  of  the  calculation: 

1 1  -  min  1 ,  ^  p  (0\  =  0m  |  A )  j 


•  )• 


17 


dominates  the  interval  error  probability  for  SNRs  in  the  vicinity  of  the  threshold  region.  These 
interval  error  probabilities  are  sample  support  and  SNR  dependent  in  general  and  essentially  weight 
the  transition  of  the  MSE  from  the  asymptotic  in  low  SNR  region  (no  information  region)  to  the 
asymptotic  in  high  SNR  region  (described  by  VD  approximation) ;  thus,  producing  the  threshold  or 
ambiguity  region  of  the  MSE  curve  as  SNR  is  varied.  An  estimate  of  the  estimation  threshold  SNR 
is  easily  obtained  from  the  predicted  MSE  curve. 


4.1.2  Multiple  Sources,  K  >  1 


Now  assume  arbitrary  K  >  1;  in  addition  assume  that  these  K  signals  are  separated  by 
at  least  a  beamwidth  (consideration  of  closely  spaced  sources  shall  be  given  in  Section  (i).  The 
extension  of  MIE  to  multiple  sources  is  accomplished  by  expanding  the  NIE  set  to  include  all  local 
neighborhoods  of  the  K  peaks  in  the  ambiguity  function  due  to  the  K  sources  present.  As  in  the 
previous  section,  the  large  sample  local  error  MSE  approximation  obtained  via  cryB{9k)  will  be  used 
to  describe  the  NIE  component  contribution  to  the  over  MSE  of  the  k-th  source  parameter  estimate. 


Let  all  local  maxima  of  the  ambiguity  function  Weapon  (0)  within  the  signal  parameter  domain 
of  interest  when  evaluated  at  K  large  SNRs  (large  enough  such  that10  9k  =  9k ,  k  =  1,2,...,  AT) 
be  given  by  the  finite  set  M  =  {  9  \  9\,  9^, • . . ,  Ok,  0/c+i,  •  •  • ,  9k+m- l },  where  9k  for  k  =  1,2,...,  K 
represent  the  peak  locations  due  to  the  K  sources,  and  9k  for  k  =  K  +  1,  A’ ’  +  2, . . . ,  A'  +  M  —  1 
represent  all  other  local  maxima.  Local  errors  for  the  k-th  source  will  originate  from  the  local 
neighborhood  of  the  9k  peak  of  the  ambiguity  function.  Thus,  signal  parameter  estimates  falling 
within  the  local  neighborhood  of  the  ambiguity  function  about  the  parameter  value  9k  will  be 
classified  as  NIE  events.  Global  errors  for  the  k-th  source  can  originate  from  all  nonsource  intervals, 
i.e.,  intervals  centered  about  local  maxima  9k  for  k  >  K.  Intersource  errors,  for  example,  error 
contributions  to  the  MSE  of  the  parameter  for  source  k  =  k\  due  to  source  k  =  ^  k\,  are  very 

unlikely  to  occur  because  the  sources  are  assumed  separated  by  at  least  a  beamwidth.* 11  Consider 
the  k-th  source  with  true  parameter  value  9  —  9k-  The  total  MSE  for  this  Capon  parameter 
estimate  can  be  approximated  by 


£|(&-0*)  \a\*  1-  Y.  p(y0k  =  Om\A) 

^  '  rn=K  + 1 


K+M—l 

+  Y  p(dk  =  0m\A)(em-ek)2. 

m=K+\ 


GVB^k) 


(18) 


The  interval  error  probability  p  y9k  =  9m |  Aj  represents  the  likelihood  of  the  Capon  search  algo¬ 
rithm  choosing  the  false  peak  located  at  9  —  9m  as  an  estimate  for  9k,  when  the  K  true  signals  are 
located  at  parameter  values  9  =  9k,  k  =  1,2,...,  K. 


Use  of  (18)  or  (17)  requires  calculation  of  the  interval  error  probabilities.  Exact  calculation 

10Such  SNRs  will  exist  provided  that  no  array  response  mismatch  is  present,  i.e.,  provided  that  the  array  responses 
used  to  compute  il>caPon(Q)  match  the  K  array  responses  existing  in  the  true  data  covariance  R  for  Ok ,  /c  =  1,2 . K. 

1 1  Intersource  errors  are  not  accounted  for  by  the  large  sample  MSE  predictions  of  VB,  nor  the  adaptation  of  MIE 
as  presented  herein,  and  does  not  apply  to  the  closely  spaced  sources  case.  However,  the  probability  of  resolution 
provided  herein  should  prove  instrumental  for  such  future  extensions. 
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is  very  difficult.  Consequently,  the  Union  Bound  (UB),  which  is  widely  used  in  Digital  Communi¬ 
cations  to  calculate  A/-ary  detection  error  probabilities  for  the  determination  of  bit /symbol  error 
rates  [5],  is  employed  to  simplify  the  analysis  while  maintaining  an  acceptable  degree  of  accuracy. 
Define  the  set  of  M  positive  integers  A4  =  {fc,  K  +  1,  K  +  2, . . . ,  K  +  M  —  1}.  As  in  [67,  2],  note 
that  the  interval  error  probability  for  the  A/-ary  hypothesis  testing  problem  is  given  by 


/ 


=  l-p(6kJ:0m\A) 


1  -  Pr  < 


U 

n  6  Afk 
n  ^  m 


[P Capon  (fin)  >  Pcapon  (^m)|  > 


where  the  II B  leads  to 


Pr 


u 

n  e 
n  ^  m 


[Pcapon  fin)  ^  Pcapon  fim)  I  *4] 


<  y  p'ipc  apon  (fin  )  ^  Pcapon  (^rn)l  *4] 

n  e  A fk 
n  /  m 


with  the  upper  bound  determined  by  pairwise  error  probabilities.  Methods  exist  to  tighten  this 
bound  by  expurgating  redundant  terms  from  the  sum  [62,  10].  A  simple  approach  that  improves 
the  UB  approximation  in  all  cases  considered  herein  is  to  approximate  the  sum  with  the  dominant 
term  [5].  This  leads  to  the  approximation 


^  ^  Pr  [  Pcapon  fin)  >  Pcaponfim)\^\  —  Pr  [  Pcapon  fik)  ^  Pcapon  fim)  |  -A] 

neA4  (19) 

n  ±  m 

=  1  ~  Pr  [Pc apon  fim)  >  Pcapon  fik ) |  -4]  ; 


thus,  the  interval  error  probabilities  are  approximated  by  the  dominant  pairwise  error  probability 


P 


=  er] 


•4)  —  Pr  [Pc  apon  fim)  >  Pc  apart  fik)  \  *4]  . 


(20) 


This  simple  approach  to  approximating  the  interval  error  probabilities  works  remarkably  well  for 
the  MSE  prediction  of  DOA  estimates  of  planewaves.12 


To  summarize,  the  goal  is  to  approximate  the  MSE  of  the  Capon  signal  parameter  estimate 
for  the  A*- tli  source  via  us< 

These  interval  error  probabilities  will  be  approximated  via  the  dominant  term  of  the  UB  sum,  as 
in  (20).  It  is  found  that  this  modified  UB  approximation  is  remarkably  accurate  in  the  vicinity  of 

1  “Initial  results  [30]  demonstrate  very  good  prediction  of  localization  performance  of  adaptive  matched  field  pro¬ 
cessing  of  underwater  acoustic  field  data  (often  a  high  multipath  environment).  Although  these  initial  results  are 
encouraging,  better  approximation  of  the  exact  interval  error  probabilities  should  be  considered  an  open  problem. 
The  simple  approach  taken  herein  is  known  to  break  down  in  the  multidimensional  parameter  case,  z.e.,  when  each 
signal  is  described  by  a  set  of  signal  parameters  Ok  =  [0/k,i,0fc.2?  •  •  •  i@k.A]T ,  although  the  idea  of  projected  ambiguity 
function  has  shown  promise  [67]  in  repairing  this  breakdown. 


3 


of  (18)  that  requires  calculation  of  error  probabilities  p  ( 0^  =  9m  ,4 ) 
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the  estimation  threshold  SNR.  but  tends  to  over  predict  the  MSE  in  the  no  information  region. 
Thus,  the  minimum  of  (18)  and  the  worse  case  MSE  obtained  with  an  estimate  Ok  that  is  uniformly 
distributed  over  the  parameter  search  space  will  be  chosen  as  the  MSE  prediction . 

Regarding  the  Bartlett  algorithm,  simply  note  that  the  sets  of  local  maxima  Ai\  =  {  6  |  0\ ,  62, 

. . . ,  0][j }  and  M  =  {0  \  0\,  62, . . . ,  0/c,  0/c+i, . . .  ,0/c+M-i}  will  be  based  on  ambiguity  function 
-0 Bartlett{0 )•  The  single  source  (K  =  1)  and  the  multiple  source  (K  >  1)  based  MSE  for  the 
Bartlett  DOA  estimates  are  given,  respectively,  by  the  following: 


E 


M 

A )  (1 9m-0i )2 

m= 2 


•  a 


H 


nW  •) 


(21) 


4)2 

+  Y, 


K+M  —  l 


K+M—l 

Y  p(0k  =  6m\A) 

m=K+\ 


(9m  -  4)2  • 


a  Tin  (9  k-) 


(22) 


The  Bartlett  interval  error  probabilities  are  approximated  by  the  following  pairwise  error  probabil¬ 
ity13 


p(j)k  =  8m  |  -4)  —  Pr  [ Pb artlett (0m)  >  PBartlett(0k)\  A\  .  (23) 

The  next  section  summarizes  the  algorithm  for  computing  the  Capon  and  Bartlett  pairwise  error 
probabilities. 


4.2  Capon  and  Bartlett  Pairwise  Error  Probabilities 


The  Capon  statistic  evaluated  at  two  test  points  is  given  by 


Pcaponifia)  — 


L-N  +  l  v"(0o)R-i  v(0a) 


Pcapon  i^b)  — 


L-N+l  vH(0b)R-lv{6b) 


(24) 


The  desired  pairwise  error  probabilities  are  of  the  form 


PeCap™(9a\9b)  =  Pr  [Pcapon(0a)  >  Pcapon(9b)\ A]  .  (25) 

1,1  Although  no  mismatch  has  been  assumed  thus  far,  when  K  >  1  the  signal  sidelobe  interference  of  multiple  signals 
introduces  mismatch  effects  (biases)  with  the  Bartlett  algorithm.  The  presence  of  these  effects  is  already  accounted 
for  in  the  HN  asymptotic  local  error  MSE  term,  i.e the  NIE  term,  but  the  IE  term  in  theory  also  requires  that  the 
error  probabilities  be  computed  with  respect  to  the  asymptotic  maxima  The  details  of  extending  analysis  to  the 
mismatch  case  is  discussed  further  in  Section  8. 
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The  following  probabilities  are  based  on  the  assumption  that  the  array  spatial  snapshots  are  in¬ 
dependent,  identically  distributed  zero  mean  complex  Gaussian  with  covariance  R.  Define  the 
following  function 


F(x,  J )  = 


(1  +  x)2J~l 


J- 1 


E 


2  J  —  1 
k  +  J 


(26) 


where  F(:r,  J )  is  the  cdf  for  a  special  case  of  the  complex  central  F  statistic;  namely,  the  case  when 
both  degrees  of  freedom  are  equal  (see  Appendix  C).  The  algorithm  for  computing  the  pairwise 
error  probabilities  for  the  Capon  estimator  is  as  follows: 


1.  Choose  the  desired  covariance  parameter  R.  Although  the  resulting  probabilities  hold  for  any 
hermitian  positive  definite  (hpd)  R,  the  following  form  satisfies  the  conditioning  event  A: 


R  =  KN  +  '£*lk-v(6k)vH(0k). 


(27) 


A:=  1 


2.  Define  parameters  P+(F).P_(F)  and  Iab{F )  such  that: 
A  vH(9b) R-lv(0b)  ±  FwIi(ea)R-1v(ea) 


P±(F )  = 

Define  the  ratio  l\(F)  and  sign  variable  S\(F): 


and  /n6(F)  =  F|vw(0fl)R^v(M  •  (28) 


P-(F)  +  ^?(F)  -  «F)  : 

l\(F)  = - -  ...  -  -==;  S\(F)  =  sign 


P-(F)~  v/p2(F)-7a6(F) 


P-(F)  -  \J P‘i{F)  -  4fr(F) 

3.  The  desired  exact  pairwise  error  probability  for  the  Capon  algorithm  is  given  by 


(29) 


Peap(m{Oa \0b)  =  0.5  •  {1  +  5a(1)}  -  SA(1)  •  F[/a(1),  L-N  +  2], 


(30) 


These  pairwise  error  probabilities  conveniently  consist  of  finite  sums  involving  no  numerical  inte¬ 
gration. 


The  Bartlett  statistic  evaluated  at  two  test  points  is  given  by 

PBartlett(Oa)  =  J  ■  VW(0a)Rv(0a),  PBartlett{Ob)  =  \  •  VW(06)Rv(06).  (31) 

The  desired  pairwise  error  probability  is  given  by  the  following 

PeBartlett(0a\6b)  =  Pr  [PuarUettiOa)  >  P liar  Betti®  b)\  -A]  .  (32) 

The  algorithm  for  calculation  of  the  Bartlett  error  probability  is  exactly  the  same  as  the  Capon 
with  minor  parameter  swaps.  In  particular,  the  parameter  definitions  of  step  2  should  replace  R  1 
with  R.  6a  with  and  Ob  with  9a,  and  in  step  3  L  —  N  +  2  should  be  replaced  with  L  in  (30).  The 
proof  of  these  probabilities  is  given  in  Appendix  A. 
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4.2.1  A  Simple  Large  Sample  Approximation  to  Pairwise  Error  Probabilities 

While  the  finite  sum  of  (30)  is  simple  to  compute  for  small  values  of  L  —  N  +  2,  it  can  become 
problematic  for  large  values  due  to  the  binomial  coefficient.  The  ratio  /a(1)  approaches  zero  as 
the  SNR  increases.  Thus,  the  pairwise  error  probabilities  of  the  Capon  algorithm  are  given  by 
the  tails  of  the  complex  F-statistic  F/,_jv+2,l-tv+2  that  are  equivalently  represented  by  the  tails 
of  an  associated  complex  beta  statistic  /?/,_/v+2,L-/v+2  (see  Proposition  C.5  of  Appendix  C).  In 

particular,  let  J  =  L  —  N  +  2,  and  note  that14  Fjj  =  —  1  +  l//3jj.  Thus,  it  follows  that  the 
probability  in  the  tails  can  be  derived  from  the  relations 

Pr  (Fjj  <x)  =  Pt  (-1  +  1/Pjj  <x)  =  Pr  6jj  >  •  (33) 

Note  that  the  pdf  of  this  beta  random  variable  has  perfect  symmetry  in  [0, 1]: 

P0j,j(0)  ocI/JU-/?)]-7-* 1.  (34) 

The  first  and  second  moments  are  easily  shown  to  be 

*«<•»> -j.  EPM  =  mTTy  <:!5> 

Due  to  the  perfect  symmetry  of  the  pdf  of  this  beta  statistic,  it  is  reasonable  to  approximate  the 
beta  tails  with  those  of  a  Gaussian  via1'’ 


1  1  \ 


Pr  (Fjj  <  x)  ~  Q 


1  + X  2 

I  ~J  + 1  r 

Vy  2(2  J  +  1)  _  4/ 


(36) 


where  Q(  )  is  related  to  the  complementary  error  function.  Since  this  error  function  is  based  on 
the  real  Gaussian  distribution  with  zero  mean  and  variance  of  one  half,  the  m  standard  deviation 
point  in  the  tail  is  given  by  satisfying  the  equation 


1  1 
1  +  x  ~  2 


J  +  1 
2(2.7+  1) 


1 

4 


m 

71' 


(37) 


Solving  for  x  one  obtains  the  desired  threshold  for  the  m  standard  deviation  tail  point: 


a(m,  J)  ~ 


1  - 


rn 


V'2JTT 


1  + 


rn 


V2IT\. 


(38) 


14 II  random  variable  A  has  the  same  pdf  as  random  variable  B.  then  they  are  said  to  be  equal  in  distribution,  and 
this  is  denoted  by  A  =  B. 

1  ’The  complementary  distribution  function  of  the  standard  Gaussian  is  denoted  by  2(  ).  It  is  related  to  the 
well- tabulated  error  function  erf(-)  and  complementary  error  function  erfc(-);  namely,  Q(x)  =  ^er^(  (^)  = 

\  |^1  —  erl  (v^f)]  =  ^57  e~*  /2dt-  See  [1,  5]  for  more  details. 
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This  tail  threshold  will  be  of  use  in  Section  5.  Equation  (36)  leads  to  the  following  large  sample 
approximation  of  the  Capon  algorithm  pairwise  error  probabilities: 


pCapcm 


(0a|06)~O.5-{l+SA(l)}-SA(l)-e 


i _ i  \ 

1  +  ^a(I)  2 


(39) 


L-N  +  3  1 


\\2[2{L-  N  +  2)  +  l]  4/ 


Equation  (36)  similarly  leads  to  the  following  large  sample  approximation  of  the  Bartlett  algorithm 
pairwise  error  probabilities: 


!  1  \ 
1+/a(1)  2 


pBartletti9a\9b)  „  Q  5  .  {1  +  SA(1)}  _  $A(1)  .  Q 


(40) 


L  +  1  1 


\  2{2L  +  1)  4  J 


where  for  the  Bartlett  algorithm  J  —  L  and  the  aforementioned  parameter  swaps  have  been  as¬ 
sumed.  These  approximations  work  well  for  large  values  of  the  degrees  of  freedom  parameter  (good 
for  J  >  40  approximately),  and  can  be  valuable  for  system  design  and  analysis,  since  the  tail 
probabilities  of  the  Gaussian  are  well  tabulated.  Indeed,  this  approximation  proves  quite  useful  for 
predicting  the  Capon  (and  Bartlett)  threshold  SNR  directly  in  closed  form. 


5.  A  METHOD  OF  DIRECT  THRESHOLD  SNR  APPROXIMATION:  A 
SINGLE  PLANEWAVE  SIGNAL  IN  WHITE  NOISE 


The  method  of  threshold  SNR  prediction  presented  in  Sections  4.1  4.2  requires  an  initial 
prediction  of  the  MSE  curve,  from  which  a  subsequent  estimate  of  the  threshold  SNR  can  be  easily 
inferred.  The  predicted  MSE  curve  is  useful,  as  it  extends  the  measure  of  parameter  accuracy  of 
the  Capon  (and  Bartlett)  algorithm  well  into  the  threshold  region.  It  is  desired,  however,  in  this 
section  to  consider  direct  calculation  of  the  threshold  SNR  that  does  not  require  an  initial  MSE 
prediction.  This  will  be  done  for  the  canonical  case  of  a  single  planewave  signal  in  white  noise. 

Steinhardt  and  Bretherton  [48]  derived  a  closed  form  expression  for  the  threshold  SNR  of 
ML  estimation  of  the  frequency  of  a  sinusoid  in  AWGN.  The  simple  analytic  forms  of  the  outlier 
probabilities  and  the  asymptotic  MSE  expression  (given  by  the  CRB)  made  such  a  calculation 
tractable.  Since  such  niceties  are  not  available  for  the  Capon  algorithm,  the  approach  taken  herein 
will  build  upon  the  premise  that  the  MSE  consists  solely  of  local  error  contributions  when  the 
search  algorithm  obtains  its  maximum  from  the  local  interval  containing  the  unique  global  maximum 
of  the  underlying  ambiguity  function  with  near  probability  1.  The  threshold  SNR,  or  equivalently 
the  smallest  SNR  at  which  this  local  error  domination  occurs,  will  be  approximated  by  using 
the  pairwise  error  probability  associated  with  the  most  significant  contender  to  the  true  global 
maximum;  namely,  the  location  of  the  largest  false  peak  of  the  ambiguity  function. 

Let  the  location  of  the  global  maximum  of  the  ambiguity  function  ipCapon(Q)  be  given  by 
Oqm  =  0\ ,  and  let  the  location  of  the  largest  competing  local  maximum  (location  of  the  highest 
sidelobe  level)  be  given  by  9  is  =  Q‘i-  The  true  data  covariance  is  clearly  given  by  R  =  ctH  + 
(t‘sv(0 \)vn (0\).  For  simplicity  of  notation,  define  the  following  variables 

x  =  a%/a~.  v\  =  ||v(#i)||2,  v2  =  ||v(02)||2,  V12  =  |vh(02)v(0i)|2 ,  C  =  l/a2.  (41) 

Using  the  matrix  inversion  lemma,  it  is  straightforward  to  show  that  the  parameters  of  step  2  of 
Section  4.2  can  be  written  in  the  form 


p±(1)  =  5 


c 

i 

x  H - 


±<>2T 


C^12 

X  + 


/li(l)  = 


C2't’l2 

(1  +  XV\)2’ 


(42) 


The  ratio  l\(  1)  is  clearly  SNR  dependent  and  defines  the  argument  of  the  F  distribution  in  (30) 
and  the  complementary  error  function  in  (39)  that  determines  the  probability  of  an  outlier.  Thus, 
the  smallest  SNR  at  which  local  error  domination  occurs  can  be  inferred  by  finding  the  value  of 
l\(l)  necessary  for  “near  enough”  to  zero  error  probability.  The  needed  proximity  to  zero  error 
probability  is  in  general  a  fairly  complicated  function  of  (i)  the  ambiguity  function  and  its  several 
derivatives  evaluated  at  the  true  source  location,  and  (ii)  the  mathematical  form  of  the  error 
probability  itself.  Empirical  analysis  of  uniform  linear  arrays  (ULAs)  indicates,  however,  that  for 
the  large  sample  case  the  error  probability  when  evaluated  at  the  threshold  SNB  is  approximately 
that,  obtained  about,  four  standard  deviations  out,  in  the  tails  of  the  Gaussian  approximation  (39). 
Denote  the  threshold  value  of  the  m  standard  deviation  tail  by  a  that  is  approximately  given  by 
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(38).  Define  the  variables 


1  +  a 
1  —  a 


1  2 


>  = 


A  (,  +  1  2  A  Vn  2  A  ,  2 


j  •  Ctf>  — 


( ' 1 1  ’2 


,  4  =  !-4’  and  llv(^)ll2  =  7,  for  k  —  1.2  (43) 


where  ci  is  the  geometric  cosine  [18]  and  it  is  assumed  that  the  scanning  vectors  are  of  equal  norm.  It 
can  be  shown  [with  (42)  and  elementary  algebra]  that  satisfying  the  constraint  l\{l)  =  o(4,  L—N+2) 
requires  obtaining  the  correct  real  root  of  the  following  quartic  equation  in  the  SNR  variable  x: 


x4  +  ax3  +  bx 2  +  cx  +  d  =  0 

where  the  coefficients  of  the  quartic  are  given  by  the  following  equations: 


(44) 


a  =  27s£  +  -(1-V>) 

,  7 

— ^ - h  72si  —  2 ci  +  2(1  —  V;)  (— rj  +  1  +  s'i 


b  = 


72(^-l) 


c  = 


lit  -  f) 


+  2 


0 
‘2  > 


(?+1+s*) 


-  (!+*«)  +l4 


-  2  if> 

2  /  1  \  1 1 

4  b  +  -  +  - 

. 

L  V  7/  7  J 

(45) 


d  =  +  4  +  2*0(1  -  if>). 

One  can  argue  that  at  least  one  of  the  roots  is  guaranteed  to  be  real  and  non-negative,  yielding 
the  desired  threshold  SNR  by  the  intermediate-value  theorem  of  elementary  calculus  when  applied 
to  the  continuous  F  distribution  of  (30),  or  its  continuous  large  sample  approximation  (39).  The 
closed  form  solution  to  the  quartic  equation  is  known  [56,  22].  The  common  approach  is  to  convert 
the  quartic  to  a  subsidiary  or  depressed  cubic: 


y3  +  py2  +  qy  +  r  =  0; 

where  p  =  6;  q  =  ac  —  4 d;  r  =  a2 cl  +  c2  —  4 bd. 


(46) 


The  known  solution  for  the  roots  of  a  cubic  equation  can  now  be  applied,  from  which  the  solution 
of  the  original  quartic  will  then  follow.  Define  the  variables 


u  —  q  —  p2/ 3;  v  =  r  —  pq/ 3  +  2(p/3)3;  D  =  4(u/3)3  +  v2. 

One  of  the  roots  of  the  subsidiary  cubic16  is  given  by  the  following  equation 

(u/3)[2/(s/D  +  i>)](1/3>  ~[(y/D  +  v)/2]W*)  -  p/3,  D  >  0; 


(47) 


Y  =  < 


2  (— ?i/3)  x  cos<  -  c 


‘OS 


—v 


2(— u/3)(3/2)_ 


(48) 


-  pi 3,  D  <  0 


16There  are  obviously  three  roots  to  the  subsidiary  cubic,  and  in  theory  it  doesn’t  matter  which  solution  is  used  to 
obtain  the  final  solutions  to  the  quartic  equation.  Thus,  the  root  possessing  the  simplest  analytic  form  was  chosen 
for  this  presentation. 


where  the  form  of  the  solution  depends  on  the  sign  of  the  discriminant  D.  Defining  the  variables 


9  =  (a/2)  -  V(a/2)2  -b-Y;  h  =  -(Y/2)  -  y/WJW  Z  d,  (49) 

the  desired  real  11011-negative  root  of  the  original  quartic  that  yields  the  desired  threshold  SNR  of 
the  Capon  algorithm  is  given  by 


SNR$P  (C2,7 ,L,N)  ~  —(g/2)  +  ^(g/2)*  -  h. 


(50) 


Applying  the  same  approach  to  the  Bartlett  algorithm  results  in  the  quadratic  equation 
Ax2  -1-  Bx  +  C  =  0,  whose  solution  yields  the  desired  threshold  SNR.  The  coefficients  are  easily 
shown  to  be 


i  =  74(£-l)4  B  =  -A^s%  C  =  — 472s|.  (51) 

The  desired  threshold  SNR  for  the  Bartlett  algorithm  uses  tail  threshold  a(4.  L)  and  is  given  by 
the  quadratic  formula: 


CMP  Bartlett  (  2  ^  r  \r\  ^  ~ 

1  +  x/1  +  *}(£-!)_ 

SNKr//  ^,7 

While  the  expression  for  the  Capon  threshold  SNR  could  benefit  from  further  asymptotic  simplifi¬ 
cation,  the  Bartlett  threshold  SNR  expression  is  rather  intuitive.  Note  that  the  Bartlett  threshold 
SNR  is  inversely  proportional  to  the  steering  vector  norm  7.  This  makes  sense  because  the  norm 
directly  impacts  the  element  level  SNR;  in  fact,  one  could  assume  unit  norm  and  absorb  7  into  the 
element  level  SNR.  Larger  7  leads  to  high  element  level  SNRs,  and  thus  a  lower  threshold  value 
can  be  tolerated.  Smaller  7  leads  to  low  element  level  SNRs.  and  thus  a  higher  threshold  value  is 
needed.  It  is  also  reasonable  to  interpret  7  as  a  measure  of  the  effective  array  gain.  Note  further 
that  the  Bartlett  threshold  SNR  is  also  inversely  proportional  to  the  geometric  sine  variable  s2. 
If  the  steering  vector  for  the  largest  sidelobe  ambiguity  direction  is  strongly  correlated  with  the 
steering  vector  of  the  true  signal  direction,  then  the  competing  sidelobe  level  will  be  high  and  s2 
will  be  small.  As  s2  becomes  small,  the  Bartlett  threshold  SNR  expression  predicts  that  a  much 
higher  threshold  SNR  value  is  necessary  to  overcome  the  influence  of  the  competing  sidelobe. 

The  utility  and  accuracy  of  these  results  is  demonstrated  in  Section  9. 
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6.  A  MEASURE  OF  THE  PROBABILITY  OF  RESOLUTION  FOR  THE 
CAPON  AND  BARTLETT  ALGORITHMS 

Empirical  analyses  (presented  in  Section  9)  corroborate  the  success  of  the  adaptation  of  MIE 
to  the  Capon  and  Bartlett  algorithms  presented  in  Section  4  for  well  separated  sources.  When 
sources  are  not  well  separated,  but  are  closely  spaced,  then  empirical  based  MSE  curves  tend  to 
depart  in  form  from  the  typical  composite  curve  illustrated  in  Figure  1.  MIE  as  described  in  Section 
4  will  not  predict  these  observed  nontypical  MSE  curves,  because  it  does  not  account  for  the  impact 
of  unresolved  sources  on  the  ambiguity  function  (at  low  enough  SNRs  unresolved  sources  lead  to  a 
single  peak  of  the  ambiguity  function),  and  it  does  not  account  for  intersource  error  contributions 
to  the  MSE,  that  can  occur  even  after  signals  are  well  resolved. 

A  useful  measure  of  the  probability  of  resolution,  however,  can  be  defined  that  provides  good 
prediction  of  the  SNR  at  which  sources  can  be  resolved  by  the  Capon  algorithm  (and  Bartlett, 
although  the  usual  Fourier/Rayleigh  limit  suffices  for  this  conventional  approach),  that  in  turn 
serves  as  a  lower  bound  on  the  threshold  SNR  [25].  Appendix  A  provides  the  derivation  of  the 
exact  cdf  and  pdf  for  the  following  ratio  of  power  spectral  estimates: 

r  ,n  n  ,  _  PcaponWa)  _  v"  (#fe)R~  *  v(^)  .  (ft  A  PBartlettjOa)  _  v"  (0„  )Rv{0a )  .  . 

Pcaponi^b)  V11  (0a)f{~]  v(0a)  ^  P Bartlett(^b)  R.v(0fc) 

The  derived  cdf  and  pdf  hold  for  any  arbitrary  hpd  data  covariance  R.  For  a  two  closely  spaced 
sources  scenario  of  the  form 

R  =  R.v  +  a2Sav(do)vH(0o)  +  <46v(0 o  +  S0)vh{9  „  +  60),  (54) 

define  parameter  Omp  as  the  parameter  value  for  the  source  with  the  smallest  power  out  of  the 
ambiguity  function,  i.e., 


0 M p  =  arg  min  tpcaP<m(0) •  (55) 

0o.eo+66 

A  two-point  measure  of  the  Capon  probability  of  resolution  can  be  defined  as 

PreTm{6O,0O  +  S0)  =  Pr 
=  Pr 

where  it  is  assumed  that  0  <  p  <  1  is  a  parameter  related  to  confidence  (probabilities  actually 
hold  for  all  non-negative  p).  The  parameter  p  essentially  defines  the  desired  “dip”  in  the  Capon 
power  spectrum  between  two  closely  spaced  sources.  For  example,  if  and  p  =  0.5,  then 

PresPon(0 0,00  +  SO)  is  the  probability  that  the  dip  in  Pcap<m(0)  midway  between  these  two  sources 
is  at  least  3  dB  less  than  the  peak  of  either  source.  Similar  measures  of  resolution  have  been 
proposed  in  [18,  25,  50,  72,  61,  47].  The  algorithm  for  computing  the  Capon  two-point  probability 
of  resolution  is  as  follows: 


Pcapon  ^0  T  2  j  —  P  '  P Capon M p) 
Pa  +  — ,  Omp^\  <  p 


(56) 
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1.  Choose  the  desired  covariance  parameter  R  as  in  (54)  and  parameters  in  equations  (28)  (2!)) 

with  F  =  p  to  obtain  P±(p),  Iab{p )<  (p),  and  /a(p)  with  0a  =  0()  +  50/2  and  6^  =  6\ip. 

2.  The  desired  two-point  probability  of  resolution  is  given  by 


This  probability  of  resolution  in  the  large  sample  case  can  be  well  approximated  as 


Pre7m{e 0. 0O  +  66)  a  0.5  •  {1  -  SA(p)}  +  Sx{p)  ■  Q 


1 


1 


1  +  Mp)  2 


L-N  +  3 


\J  2[2(L  -  TV  +  2)  +  1]  A) 


The  Bartlett  probability  of  resolution  is  obtained  from  a  similar  formulation: 


PrBertlett(0o,  0o  +  69)  =  Pr 
=  Pr 


P Bartlett  (00  +  <  P  '  P, Bartlett  {6 Mp) 

h(t>o+d4-0Mp)  <p 


(57) 


(58) 


(59) 


The  algorithm  for  computing  the  Bartlett  probability  of  resolution  is  very  similar  to  that  for 
computing  the  Capon  probability  of  resolution  with  small  parameter  changes: 


1.  Choose  the  desired  covariance  parameter  R  as  in  (54)  and  parameters  in  equations  (28)  (29) 
with  F  =  p  to  obtain  P±(p ),  Iab(p )•  S\ (p),  and  l\{p)  with  0a  =  6mp  and  0^  =  0o  +  50/2  and 
R“ 1  replaced  by  R . 

2.  The  desired  two-point  probability  of  resolution  is  given  bv 


This  probability  of  resolution  in  the  large  sample  case  can  be  well  approximated  as 


i 


PrBe?tlett{6 o.  0O  +  450)  ~  0.5  •  { 1  -  Sx(p)}  +  Sx(p)  ■  Q 


1  \ 


VN 


i  +  i\{p)  2 


L  T  1  1 


2(2L  +  1)  4/ 


(00) 


(61) 


The  numerical  results  of  Section  9  will  demonstrate  the  utility  of  these  measures  for  predicting 
the  SNRs  at  which  closely  spaced  sources  become  resolvable  by  the  Capon  or  Bartlett  method.  This 
resolution  SNR  also  provides  a  lower  bound  for  the  threshold  SNR  above  which  the  large  sample 
MSE  performance  predictions  of  VB  and  others  apply  when  closely  spaced  sources  are  present  . 
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7.  DIRECT  CLOSED  FORM  CALCULATION  OF  THE  CAPON 
RESOLUTION/DETECTION  SNR:  TWO  EQUAL  POWER  SIGNALS  IN 

WHITE  NOISE 


The  goal  of  this  section  is  to  provide  a  method  for  computing  the  Capon  resolution  SNR 
directly  that  does  not  require  an  initial  calculation  of  the  probability  of  resolution  curves.  This  can 
be  accomplished  in  a  manner  similar  to  that  used  for  direct  approximation  of  the  threshold  SNR 
described  in  Section  5. 

Let  the  desired  probability  of  resolution  be  given  by  Pr/,  e.g.,  if  a  90%  likelihood  of  a  3 
(IB  dip  (p  =  0.5)  occurring  between  two  spectral  peaks  is  desired,  then  Pd  =  0.90.  This  desired 
resolution  probability  implies  a  necessary  threshold  value  for  the  ratio  l\(p).  For  small  sample 
support  L,  or  in  particular  small  values  of  L  —  N  +  2,  one  must  numerically  invert  the  F-distribution 
given  in  (57)  to  obtained  the  needed  threshold  d.  Such  inversions  are  commonplace  in  adaptive 
detection  analysis,  as  they  provide  the  required  threshold  values  for  detection  statistics  to  maintain 
a  specified  probability  of  false  alarm  [27,  26,  49].  For  large  values  of  L  —  N  +  2  (approximately 
greater  than  40),  however,  a  simple  inversion  is  possible  using  the  error  function.  Recall  that  the 
probability  is  determined  by  the  tails  of  the  complex  F  distribution,  i.e.,  a  calculation  of  the  form 
Pr  [FJyJ  <  d(c;  J)]  =  e.  Using  the  large  sample  approximation  of  the  F  distribution  discussed  in 
Section  5  and  presented  for  the  probability  of  resolution  in  equation  (58)  of  Section  6,  it  is  easy  to 
show  that  the  desired  tail  threshold  is  well  approximated  by  inverting  the  error  function  to  obtain: 


a(c,  ,J) 


T±_L  _  I  j  erfc"1  (2e) 
2.7+1  2  v  ' 


-l 


1. 


(G2) 


Define  the  following  variables 

e  =  0.5  •  [1  -  S\(p)]  +  S\(p)  ■  Pa,  a  =  a(e;L-N  +  2),  £  = 


d  +  1 
d  —  1 


(63) 


and  let  the  data  covariance  for  two  equal  power  closely  spaced  sources  in  white  noise  be  defined  as 

R  =  a2 1  +  <7|v(0i)vw(0i )  +  <t|v(^2)v//(6>2)  (64) 


where  clearly  the  signal  directions  are  given  by  0\  =  0q  and  0>  =  do  +  SO.  Let  the  midpoint  angle 
be  denoted  by  O3  =  0q  +  SO/2,  and  define  also  the  following  variables 

x  =  <7l/(T 2,  C  =  1/<t2,vi  =  ||v(6>i)||2,  v2  =  ||v(02)||2,  1%  =  ||v(03)||2, 
vy2  =  vH(0i)v(O2),  V13  =  v" (0i)v(0:i),  v23  =  v"(02)v{0:i) 

[note  that  U12  is  defined  differently  here  than  in  (41)].  Setting  l\(p)  =  d  and  solving  for  the  SNR 
variable  indicated  by  x,  it  can  be  shown  [with  elementary  algebra]  that  the  desired  resolution  SNR 
is  obtained  as  a  root  of  the  following  quartic  equation: 

x4  +  ax3  +  bx 2  +  cx  +  d  =  0  (66) 
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where  the  coefficients  are  given  by  the  following  ratios 


a  =  B/A,  b  =  C/A,  c=D/A.  d  =  E/A 


(67) 


that  are  defined  in  terms  of  the  following  variables 


A 

B 

C 

D 

E 


lA2_-A\  +  A\A,\2p 
2 lA-B-  -  2A+B+  +  8pRe  ( A,B ,*) 

2 lA-C-  +  B2_-  (2A+C+  +  B\)+Ap 


2Re  (A/C[)  +  \B[\2 


2£B-C-  -  2B+C+  +  SpRe  (B/C’J) 
ICl-C2  +4p|C/|2. 


(68) 


These  coefficients  are  expressed  in  terms  of  further  coefficients  defined  as  follows: 

A±  =  (i  -  («,,  ±  T  „M!  ±  2Re  ( T  (69) 

V  V\V2  I  V2  V\  \  V\V2  )  V2 


r,  /  1  ,  1  \  vl  _  _>13|2  „  _>23p  ^  (yi±pvZ) 

B±  = - 1 - ( vi  ±  pv 3) - T  P - =F  P - ,  C±  = -  (70) 

1  V\  V2 )  V2  V\V2  V\V2  V\V>2  V\V2 


1 2 


Ai  =  vi3  (  1  - 


kl2|  l  ,  ^13  1^12 n  _  (  1  ,  1 

-  I  —  ^13  H - ,  Bj  —  1>13  I - 1 - 

V\V2  I  V\V2  \V\  V2 


V2 


C/  =  (71) 

V\V2  V\V2 


The  roots  of  (66)  can  be  obtained  using  the  algorithm  described  by  (46)  (49)  from  Section  5 
swapping  (a,  6,  c,  d)  with  (d ,  6,  c,  d).  The  desired  real  non-negative  root  of  the  original  quart ic  that 
yields  the  desired  resolution  SNR  of  the  Capon  algorithm  for  confidence  level  Pd  is  given  by 


SNR^T"  (Vl,  v2,  v3,  V12,  v13,  v23.  L.  N,  p,  Pd)  =  ~{g/2)  +  y/{g/2)*  -  h. 


(72) 


Although  this  expression  could  likewise  benefit  from  an  asymptotic  simplification  [a  simple  exercise 
left  to  the  ambitious],  it  is  exact  and  found  to  be  quite  useful  for  system  parameter  trade  studies 
and  design. 
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8.  SIGNAL  MODEL  MISMATCH  EFFECTS:  MSE  PREDICTION  AND 

RESOLUTION 

The  scanning  vectors  v(0)  used  to  form  the  Capon  spectra  Pcapon(6)  are  chosen  by  the  user 
and  are  therefore  subject  to  errors.  The  array  responses  of  the  signals  present  in  the  data,  as 
described  by  the  true  underlying  covariance  parameter  R,  are  likely  to  differ  from  the  scanning 
vectors  chosen  by  the  user.  This  difference  is  historically  referred  to  as  signal  model  mismatch  [18], 
and  its  presence  often  limits  achievable  performance.  Accounting  for  its  impact  is  of  particular 
interest  for  high-resolution  approaches,  since  signals  that  do  not  fall  on  the  scanning  search  manifold 
will  suffer  some  degree  of  suppression  by  construction1 '  [see  (3)].  In  this  section,  a  strategy  for 
accounting  for  mismatch  in  MSE  prediction  is  presented  for  a  single  signal  in  white  noise.  Extension 
to  colored  noise  [30]  and  multiple  signals  is  transparent. 

Let  the  true  data  covariance  be  given  by  R  =  a2 1  +  cr|dd^,  where  the  true  signal  array 
response  is  given  by  d.  yet  the  scanning  vectors  chosen  by  the  user  are  given  by  v(9).  It  will  be 
assumed  that  although  mismatch  is  present,  the  ambiguity  function  -0 Capon(6 )  =  l/v//(^)R  -M  e) 
remains  multimodal  and  has  a  unique  global  maximum  in  the  search  domain  of  interest,  the  location 
of  which  shall  be  denoted  by  6cm •  Define  the  conditioning  event  B  =  {  Global  maximum  of 
ambiguity  function  is  at  6 cm  }•  Let  the  locations  of  all  other  local  maxima  of  the  ambiguity 
function  within  the  signal  parameter  domain  of  interest  when  evaluated  at  a  sufficiently  large 
SNR18  (large  enough  that  the  global  maximum  is  evident  and  located  at  6cm)  be  given  by  the 
finite  set  A4  =  { 6  1 0i,  62,  • . . ,  $a/-i  }.  MSE  is  defined  relative  to  some  assumed  “true”  location  of 
the  signal.  Let  that  assumed  true  location  be  denoted  6j\  Note  that  the  VB  asymptotic  local  error 
MSE  predictions  [57]  hold  for  arbitrary  data  covariance  matrices  R.  Thus,  the  asymptotic  MSE 
prediction  cryB{6r)  holds  in  the  presence  of  signal  model  mismatch,  and  no  new  theory  need  be 
developed  in  the  regime  of  high  SNR  values  for  which  these  predictions  hold.  The  Capon  asymptotic 
local  error  MSE  of  the  estimate  6  in  the  presence  of  mismatch  is  given  by 

o2vb{6t)  =  E  |  (0-  0T)2}  =  [^  (^ -  6»r}] 2  +  E  j  (o  -  0Gm)2}  (73) 

where  the  asymptotic  bias  can  be  written 

E  {fl -  07  j  =  E  {fl -  Ogm }  +  (| Ogm  -  eT) .  (74) 

Regarding  the  threshold  region  of  MSE  performance,  the  lesson  learned  in  Section  4.1  is  that 
the  correct  weights  governing  the  transition  of  MSE  performance  from  the  no  information  region  to 

1 1  While  it  is  common  practice  to  diagonally  load  the  estimated  data  covariance  to  robustify  the  Capon  beamformer 
to  the  deleterious  effects  of  mismatch  errors  [19,  61,  23],  the  impact  of  such  regularization  shall  not  be  considered 
herein.  Such  an  analyses  is  currently  underway  [42].  It  is  noteworthy,  however,  to  state  herein  that  the  asymptotic 
MSE  performance  for  a  diagonally  loaded  Capon  estimator  has  been  observed  to  be  bounded  by  that  predicted 
for  the  unloaded  Capon  estimator  [57],  and  that  predicted  for  the  unloaded  Bartlett  estimator  [21].  It  is  perhaps 
straightforward  to  construct  a  proof  of  these  bounds  based  solely  on  continuity. 

lsFor  the  DO  A  estimation  problem,  the  location  of  these  maxima  is  weakly  dependent  on  SNR,  but  in  general  this 
dependence  can  be  significant  and  should  be  accounted  for  (see  [30],  for  example). 
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the  asymptotic  region  are  those  given  by  the  interval  error  probabilities  relative  to  the  local  interval 
of  the  global  maximum  of  the  ambiguity  function.  Thus,  the  total  MSE  for  this  Capon  parameter 
estimate  can  be  approximated  by 


A/-1 

i- 

m=  1 


A/-1 


<tIb(0t)  +  (#  =  *»»|  B)  (0m  -  °r)2  (75) 


m=  1 


where  the  interval  error  probabilities  are  approximated  by  the  dominant  pairwise  error  probability 
of  the  UB  sum: 


~  Pr  [PcaponfQm)  >  Pcapon (#GA/)  |  B]  . 


(76) 


Note  that  the  estimation  error  is  with  respect  to  the  assumed  true  target  angle  6t,  but  the  error 
probabilities  are  with  respect  to  the  location  of  the  true  global  maximum  of  the  ambiguity  function 
Ogm •  Equations  (73)  (76)  represent  the  first  generalization  of  MIE  MSE  prediction  to  encompass 
signal  model  mismatch.  These  equations  provide  accurate  predictions  of  mismatch  performance 
and  have  been  recently  applied  with  success  to  fairly  complex  high  multipath  environments  where 
mismatch  is  often  inevitable  [30].  A  recent  analysis  has  lead  to  a  lower  bound  on  MSE  performance 
when  signal  mismatch  is  present  within  a  Bayesian  estimation  context  [67,  70].  The  approach 
presented  in  this  section,  although  applied  to  the  Capon  estimator,  is  likewise  applicable  to  other 
nonlinear  estimation  schemes  including  the  Bartlett  algorithm  and  ML  estimation. 

Regarding  resolution  in  the  presence  of  mismatch,  the  two- point  measure  of  the  probability 
of  resolution  derived  in  Section  6  holds  for  all  choices  of  scanning  vectors  v(0)  and  hpd  covariance 
matrices  R.  Thus,  no  modification  or  new  theory  need  be  developed  to  handle  the  impact  of  signal 
model  mismatch  on  the  resolution  of  the  Capon  algorithm. 


9.  NUMERICAL  RESULTS 


9.1  Single  Broadside  Signal  in  White  Noise 

Consider  a  DOA  estimation  scenario  involving  a  single  planewave  source  in  AWGN  resulting 
in  a  set  of  signal  bearing  snapshots  x(/)  ~  CAf[0, 1  +  (t§v(9t)vh {Ot)],  /  =  1, 2, . . . ,  L,  for  an  A  =  18 
element  uniform  linear  array  (ULA)  with  slightly  less  than  lambda  over  two  element  spacings.  The 
array  has  a  3  dB  beamwidth  of  7.2  degrees,  and  the  desired  target  signal  is  arbitrarily  placed  at 
0T  =  90  degrees  (array  broadside).  The  signal  parameter  search  space  of  interest  is  defined  to  be 
0  G  [60°,  120°].  The  signal  parameter  to  be  estimated  is  clearly  the  scalar  angle  of  arrival  0  =  0j. 
The  resulting  Monte  Carlo  simulation- based  MSE  performance,  the  Cramer-Rao  Bound  (CRB), 
and  the  MIE-based  MSE  predictions  as  a  function  of  the  array  element  level  SNR  are  shown  in 
Figure  3  for  the  sample  support  cases  of  L  —  1.5 A,  2 A,  and  3 A.  The  VB  asymptotic  local  error 
MSE  prediction  is  illustrated  for  the  L  —  2A  case  for  comparison.  Note  that  the  VB  prediction  is 
only  valid  above  the  estimation  threshold  SNR.  The  MIE  MSE  prediction  is  remarkably  accurate 
well  into  the  estimation  threshold  region,  as  it  continues  with  accurate  prediction  where  the  VB 
MSE  predictions  begin  to  become  inaccurate  at  low  SNRs.  Note  that  UB  approximation  begins 
to  over-predict  the  MSE  in  the  no  information  region.  Thus,  it  is  hard  limited  to  the  variance  of 
a  uniform  distribution.  It  is  known  that  the  Capon  estimator  is  not  asymptotically  (SNR— ►  oo) 
efficient  for  fixed  L,  such  that  increasing  the  SNR  does  not  bring  its  MSE  performance  (‘loser  to 
the  CRB  [61].  This  nonefficiency  is  likewise  manifest  in  these  performance  curves.  [MATLAB  code 
is  provided  in  Appendix  D  for  the  benefit  of  the  interested  reader,  that  will  reproduce  the  MSE 
prediction  of  this  example.] 

Figure  3  also  illustrates  the  accuracy  of  the  direct  closed  form  approximation  of  the  threshold 
SNR  given  by  equation  (50).  The  threshold  SNRs  indicated  by  the  solid  vertical  lines  use  a  value 
of  a  based  on  the  tails  of  equation  (30)  (approximately  5xl0~5  error  probability),  whereas  the 
dashed  vertical  lines  use  a  value  of  a  derived  from  the  tails  of  the  large  sample  approximation 
(39)  four  standard  deviations  out.  The  threshold  SNR  based  on  (39)  becomes  more  accurate  as  L 
increases.  Figure  5  further  illustrates  the  utility  of  this  direct  calculation  by  plotting  the  Capon 
threshold  SNR  as  a  function  of  the  peak  sidelobe  level  r*20  and  the  norm  of  the  array  responses  7. 
Smaller  norms  and  larger  sidelobe  levels  result  in  higher  threshold  SNR  points  with  the  stronger 
dependence  being  011  the  norm  parameter  7.  The  second  plot  in  Figure  5  shows  the  impact  of 
increased  sample  support  on  the  resulting  threshold  SNR  point;  reduced  finite  sample  effects  result 
in  smaller  threshold  SNRs,  with  the  stronger  dependence  being  the  height  of  the  peak  sidelobe  level 

Figure  4  shows  the  analogous  MSE  predictions  and  Monte  Carlo  results  obtained  from  process¬ 
ing  the  same  data  through  the  conventional  Bartlett  algorithm.  The  MIE-based  MSE  predictions 
are  very  accurate  well  into  the  threshold  region.  Note  that  for  this  case  of  a  single  planewave  signal 
in  AWGN  that  this  conventional  beamforming  approach  results  in  an  efficient  estimate  of  the  DOA 
at  high  SNRs.  It  can  be  shown  that  the  Bartlett  algorithm  in  fact  yields  the  ML  estimate  in  this 
case  [3,  13].  The  losses  in  parameter  accuracy  seen  by  the  Capon  algorithm  relative  to  the  CRB 
in  this  case  are  adaptive  training  losses.  Clearly,  if  there  is  only  a  single  signal  present  in  a  white 
noise  environment,  then  there  is  no  need  to  adapt  (ie.,  there  is  110  interference  to  cancel).  Since  the 
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Capon  algorithm  unnecessarily  attempts  to  adapt  in  this  benign  environment,  there  is  some  loss 
in  performance.  These  losses,  however,  become  gains  in  the  presence  of  colored  noise,  or  multiple 
signals.  Note  also  from  Figure  4  the  remarkable  accuracy  of  the  direct  threshold  SNR  calculation 
of  (52).  The  utility  of  this  direct  threshold  SNR  approximation  is  further  illustrated  in  Figure  6 
where  it  is  plotted  versus  sidelobe  level,  sample  support,  and  steering  vector  norm;  in  contrast  with 
Figure  5,  the  Bartlett  algorithm  is  the  preferred  choice  in  this  case  since  both  optimal  accuracy  is 
achieved  at  high  SNRs.  and  the  threshold  SNR  values  are  lower. 


9.2  Two  Signals  in  White  Noise 

Next  consider  the  same  scenario  as  in  Section  9.1,  but  with  an  additional  source  of  equal 
power  included  in  the  environment  at  70  degrees.  The  Capon  algorithm  MSE  performance  of  both 
signal  parameter  estimates  is  illustrated  in  Figure  7,  and  the  Bartlett  algorithm  MSE  performance 
in  Figure  8.  The  MIE  predictions  remain  quite  accurate  well  into  the  threshold  region  in  both 
cases.  This  example  also  provides  a  good  illustration  of  the  benefit  of  adaptive  processing.  The  two 
signals  have  sidelobe  energy  that  is  permitted  to  interfere  with  the  mainlobe  energy  of  the  other 
signal  when  using  the  conventional  Bartlett  algorithm.  This  interference  introduces  a  bias  (a  small 
one  in  this  case)  in  DOA  estimates  at  high  SNRs  that  cannot  be  overcome  by  larger  signal  strength. 
The  source  of  this  bias  can  be  traced  to  the  locations  of  the  peaks  of  the  ambiguity  function  relative 
to  the  true  signal  locations.  This  is  illustrated  in  Figure  9,  where  close  examination  reveals  the 
difference  in  these  locations,  and  Figure  10  shows  the  MSE  curve  of  one  signal  at  high  SNRs 
leveling  out  at  the  bias.  The  adaptive  Capon  algorithm  on  the  other  hand  places  nulls  in  signal 
directions  it  is  not  steered  to.  This  adaptive  nulling  minimizes  the  impact  of  one  signal  on  the 
DOA  estimation  of  another.  The  benefit  of  this  adaptivity  is  illustrated  in  Figure  7  as  larger  signal 
strength  results  in  improved  accuracy.  It  is  noteworthy  that  the  Bartlett  algorithm  provides  better 
accuracy  at  the  low  SNRs  in  the  threshold  region,  whereas  the  Capon  algorithm  MSE  performance 
eventually  equals  and  surpass  that  of  the  Bartlett  at  high  SNRs.  While  neither  the  Bartlett  nor  the 
Capon  algorithm  are  efficient  in  this  case,  increased  signal  strength  results  in  improved  parameter 
estimation  accuracy  for  the  Capon  algorithm. 


9.3  Mismatch:  Tilted  Minimum  Redundancy  Linear  Array 

Section  8  presented  a  strategy  for  accounting  for  mismatch  in  MSE  predictions.  As  an  illus¬ 
tration,  consider  an  N  =  4  element  minimum  redundancy  linear  array  (MR LA)  [61]  that  has  much 
higher  sidelobes  (ambiguities)  than  a  fully  populated  7-element  ULA  (see  conventional  beampat- 
terns  in  Figure  11).  Mismatch  can  be  introduced  by  tilting  this  array,  yet  processing  the  data  as 
if  the  array  were  not  tilted.  The  adaptation  of  MIE  described  in  Section  8  can  account  for  this 
mismatch  as  illustrated  in  Figure  12,  where  several  tilt  angles  have  been  considered.  Note  that  the 
MSE  predictions  accurately  capture  both  the  global  errors  due  to  the  high  ambiguities  (threshold 
region),  as  well  as  the  asymptotic  bias  resulting  from  mismatch. 
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Figure  3.  Single  source  Capon  MSE  performance ,  9\  =  0r  =  90°,  L  =  1.57V,  27V,  37V. 
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Figure  4-  Single  source  Bartlett  MSE  performance ,  0 \  —  9t  —  90°,  L  =  1.57V,  27V,  37V. 
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Figure  6.  Bartlett  threshold  SNR  for  single  source,  N  =  18  element  ULA.  and  0\  =  0t  =  90°. 
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Figure  7.  Two-source  Capon  MSE  performance ,  1500  Monte  Carlo ,  0\  =  90°,  #2  —  70°,  L  =  1.5AT,  3iV. 


Figure  8.  Two-source  Bartlett  MSE  performance ,  1500  Monte  Carlo,  0\  =  90°,  82  =  70°,  L  =  1.57V,  37V. 
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Figure  9.  Bartlett  ambiguity  function  at  high  SNR  showing  bias  origin,  N  =  18  element  ULA,  L  =  1.57V 
snapshots. 
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Figure  10.  MSE  of  Bartlett  DOA  estimate  for  N  =  18  element  ULA.  showing  tail  due  to  bias. 
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Figure  12.  Mismatch  example:  MSE  for  tilted  MRLA  with  signal  at  broadside  and  L  =  3 N.  1500  Monte 
Carlo. 
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9.4  Capon  Algorithm  Probability  of  Resolution 


The  two-point  measure  of  the  Capon  probability  of  resolution  proposed  in  Section  6  allows 
one  to  predict  the  SNRs  required  for  a  dip  in  the  estimated  power  spectrum  of  a  specified  level  to 
appear  between  closely  spaced  sources  with  high  confidence.  Consider  the  N  =  18  element  IJLA  of 
the  example  in  Section  9.2,  and  choose  a  sample  support  of  L  =  3 N.  It  is  of  interest  to  consider 
resolution  performance  as  the  position  of  the  additional  source  of  equal  power  varies.  The  resulting 
probabilities  of  resolution  for  both  the  Capon  and  Bartlett  algorithm  derived  in  Section  6  are  plotted 
in  the  upper  right  of  Figure  13  as  a  function  of  the  angle  separation  for  varied  choices  of  element- 
level  SNR.  The  superiority  of  the  Capon  method  is  clearly  illustrated  as  resolving  power  increases 
with  SNR;  in  contrast,  the  conventional  Bartlett  method  cannot  exceed  the  Fourier/ Rayleigh  limit 
no  matter  how  strong  the  signals.  This  is  reflected  in  the  Bartlett  probability  of  resolution  as  all 
curves  are  so  close  they  appear  as  one,  only  achieving  good  resolution  at  separations  of  a  beamwidth 
or  greater.  Note  that  as  the  signals  become  arbitrarily  close,  the  two-point  probability  of  resolution 
must  migrate  back  to  0.5  as  the  two  signals  become  indistinguishable. 

As  the  next  example,  let  the  additional  source  of  equal  power  be  placed  at  93  degrees,  ie., 
at  less  than  half  a  beamwidth  separation.  Let  the  nominal  ULA  array  sensor  positions  z„  be 
independently  perturbed  by  zero  mean  spherically  symmetric  Gaussian  noise,  such  that  the  true 
array  sensor  positions  are  given  by  zn  +  en.  where  e7?  ~  A/^O.  I.^^a/s)-  The  lower  left  and  right 
of  Figure  13  show  a  plot  of  the  Capon  probability  of  resolution  as  a  function  of  the  element  level 
SNR  with  L  =  2 N  spatial  snapshots.  Two  forms  of  mismatch  are  considered:  (i)  deterministic, 
and  (ii)  stochastic.  The  lower  left  plot  was  generated  by  simply  using  a  single  realization  of  a 
perturbed  array  as  the  root  MSE  of  the  perturbations  is  steadily  increased  (different  perturbations 
were  used  for  each  value  of  orjiAIS).  The  lower  right  plot  represents  the  resolution  observed  as  one 
averages  over  an  ensemble  of  perturbation  realizations.  Note  that  the  location  of  the  SNR  at  which 
resolving  occurs  appears  to  be  relatively  constant  as  a  function  of  the  perturbation  level;  namely, 
the  maximal  resolution  is  achieved  at  an  element  SNR  of  approximately  15  dB  (about  40  dB  array 
SNR).  SNR  in  excess  of  this  value  is  essentially  wasteful.  As  more  mismatch  is  introduced,  the 
asymptotic  likelihood  of  resolution  decreases.  The  dip  in  Pre^pan  before  its  eventual  rise  as  SNR  is 
varied  is  simply  due  to  the  fact  that  the  underlying  Capon  ambiguity  function  has  a  single  peak 
midway  between  the  two  sources  for  low  SNR  values  (sources  are  unresolved),  as  illustrated  in  the 
upper  left  of  Figure  13.  It  is  only  when  the  source  SNRs  exceed  that  necessary  for  resolution  that 
two  distinct  peaks  appear  at  the  source  locations,  with  a  dip  in  between.  Since  the  two-points 
used  for  the  probability  of  resolution  calculation  are  chosen  such  that  one  is  at  one  of  the  source 
locations  and  the  other  midway  between  the  two  sources,  a  dip  occurs  before  the  rise. 

Lastly,  consider  the  same  array  configuration  illuminated  by  two  equal  power  sources.  It  is 
desired  to  determine  the  SNR  required  for  these  sources  such  that  a  realization  of  the  Capon  spec¬ 
trum  will  have  a  dip  (of  any  level,  z.e.,  p  —  1)  between  the  two  signal  DOAs  with  a  90%  likelihood, 
i.e.,  Pr(sPon{0 o,#o  +  SO)  =  0.9.  This  resolution  SNR  can  be  determined  from  the  probability  of 
resolution  curves  obtained  from  (57).  One  can  likewise  make  use  of  the  direct  calculation  provided 
by  (72).  Figure  14  plots  this  resolution  SNR  as  a  function  of  both  angle  separation  60  and  sample 
support  L.  Clearly,  angle  separation  is  the  primary  determinant  for  this  resolution  SNR.  Increasing 
sample  support  provides  the  largest  gains  when  the  signals  are  appreciably  separated,  thus  driving 
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down  the  large  variance  of  the  adaptive  sidelobes  [26.  40.  41].  This  is  evidenced  in  Figure  14  by 
the  larger  swings  of  the  —5  dB  and  —10  dB  contours  as  L  is  increased.  As  the  signal  separation 
decreases,  these  swings  are  barely  noticeable  as  finite  sample  effects  are  less  important  because  the 
SNRs  required  for  resolution  are  so  large  that  they  dominate  the  convergence  of  the  adaptive  nulls 
of  the  Capon  processor.  Reference  [41]  explores  the  statistics  of  adaptive  nulling,  and  demonstrates 
that  while  convergence  is  impacted  by  finite  sample  effects,  it  is  also  inversely  proportional  to  signal 
strength.  Hence,  only  a  few  training  samples  of  a  very  strong  interferer  are  required  for  adequate 
null  formation,  whereas  many  samples  are  required  to  drive  down  the  quiescent  sidelobes  of  the 
adaptive  beamformer.  These  characteristics  of  adaptive  beamforming  are  well  known  and  serve  as 
catalysts  for  the  study  of  robust  beamforming  methods  [19.  23]. 


Figure  13.  Capon  and  Bartlett  probability  of  resolution  for  N  =  18  element  ULA. 
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Figure.  14 •  Direct  calculation  of  Capon  resolution  SNR ,  P^easpon  —  0.9. 


10.  CONCLUSIONS 


The  method  of  interval  errors  (MIE)  has  been  successfully  adapted  and  extended  to  the 
Capon-MVDR  algorithm  and  the  Bartlett  algorithm,  providing  remarkably  accurate  prediction  of 
the  MSE  threshold  SNRs  for  the  DOA  estimation  of  an  arbitrary  number  of  well  separated  sources. 
These  SNRs  are  predicted  via  simple  finite  sum  expressions  for  the  pairwise  error  probabilities, 
involving  no  numerical  integration,  and  circumventing  the  need  for  many  time-consuming  Monte 
Carlo  simulations.  Simple  large  sample  approximations  for  these  error  probabilities  based  on  the 
well-tabulated  complementary  error  function  were  derived.  Closed  form  approximations  of  the 
threshold  SNR  points  of  both  the  Capon  and  Bartlett  algorithms  were  derived  for  the  large  sample 
canonical  case  of  a  single  planewave  signal  in  white  noise.  A  new  two-point  measure  of  the  Capon 
algorithm  probability  of  resolution  was  proposed  that  accurately  predicts  the  SNRs  necessary  for 
mutual  source  resolvability  for  sources  of  arbitrary  closeness.  For  convenience,  a  direct  inversion  of 
these  resolution  probabilities  yielding  the  required  resolution/detection  SNR  was  provided  for  the 
canonical  case  of  two  equal  power  signals  in  white  noise.  This  direct  inversion  simplifies  system 
parameter  trade  studies.  Both  the  threshold  region  MSE  predictions  and  resolution  probabilities 
account  for  colored  noise,  finite  sample  effects,  and  signal  model  mismatch,  and  thus,  represent 
valuable  design  and  analysis  tools  for  any  system  employing  the  Capon  and/or  Bartlett  algorithm. 
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APPENDIX  A 

DERIVATION  OF  PAIRWISE  ERROR  PROBABILITIES  FOR  THE  CAPON 

AND  BARTLETT  ALGORITHMS 


The  analysis  presented  in  this  report  has  assumed  that  each  of  the  K  signals  of  interest 
possess  a  single  parameter  of  interest  0^.  Several  applications  involve  signals  that  possess  multiple 
parameters  each,  z.e.,  0^  =  [q^i,  &k,2, . . .  ,a^A\r  (see  [21,  24]  ).  The  following  derivation  is  valid 
for  such  a  case.  This  is  emphasized  by  the  chosen  notation. 

The  desired  pairwise  error  probabilities  are  defined  by  the  following 

P?ap0n(da \6b)  =  Pr  [P Capon a)  >  Pcapon{Ob)\  A] 

pBarUett{da  ^  =  Pr  [PBarUett(da)  >  PBartlett{Ob)\ A]  . 


The  conditioning  event  A  is  satisfied  via  the  appropriate  choice  of  the  data  covariance  R.  The 
following  derivations  are  valid,  however,  for  any  arbitrary  R:  hence,  the  conditioning  notation  is 
omitted.  The  analysis  in  [37]  requires  that  the  pdf  of  the  ratio  Fa  be  derived  where 


n  a  vH(eb)k-\(eb)  t  -  a  v"(0(,)Rv(0a) 

r  a  = - ~ -  and  F  —  - - - . 

vw(0„)R-]v(0a)  v"(0,,)Rv(06) 


(A-2) 


Both  of  these  random  variables  are  non-negative  with  probability  1,  and  hence  are  defined  over  the 
domain  of  non- positive  real  numbers.  Note  that  the  corresponding  cdfs  are  given  by 


Pr  {Pa  <  F)  =  Pr  [vH(06)R-'  v(06)  -  F  •  vw(0rt)R-1v(0(,)  <  C)] 

Pi(fa<f)  =  Pr  v"(0a)Rv(0a)  -  F  ■  v//(06)Rv(0,,)  <  0  . 

The  desired  pairwise  error  probabilities  can  be  written  in  terms  of  these  cdfs: 

PeCapo,,(0„ |0h)  =  1  -  Pr(FA  <  1)  and  PeBartlett(en\Ob)  =  1  -  Pr  (fa  <  l)  .  (A-4) 


Thus,  the  goal  of  this  section  is  to  first  derive  the  cdfs  of  Fa  and  F^  from  which  both  the  desired 
pairwise  error  probabilities  follow,  and  the  pdfs  of  these  random  variables  via  differentiation.  This 
analysis  draws  heavily  upon  the  techniques  described  in  Section  3-5  of  [49],  involving  whitening, 
ort  hogonal  rotations,  use  of  matrix  lemmas,  and  complex  Gaussian-related  statistics  (as  summarized 
in  Appendix  C).  One  such  complex  Gaussian-related  statistic  is  the  complex  central  F-statistic, 

obtained  from  the  ratio  of  two  complex  central  chi-squared  statistics,  ie.,  Xj/x%  ^  Fj.b  with  pdf 
given  by19 


PF,b(F )  = 


(J  +  B -1)1 


FJ- 


(.7-  1)!(B—  1)!(1  +  F)(J+B)’ 


F  >  0. 


(A-5) 


1JIf  random  variable  A  hai>  the  same  pdf  a.s  random  variable  B ,  then  they  are  said  to  be  equal  in  distribution,  and 
this  is  denoted  by  A  =  B. 


49 


This  statistic  will  be  useful  in  the  developments  to  follow.  The  cdfs  and  pdfs  of  Fa  and  F^  will  be 
derived  independently  in  sequence.  Let  the  N  x  L  data  matrix  and  associated  N  x  N  un-normalized 
estimated  data  covariance  from  which  the  estimated  signal  parameters  shall  be  derived  be  given, 
respectively,  by 

X  =  [x(l)|x(2)|  •  •  •  |x(L)]  and  R  =  XX".  (A-6) 


A.l  Capon  Algorithm 


Note  that  the  data  matrix  X  can  be  represented  as  a  linearly  transformed  white  data  matrix 
Z  that  is  invariant  to  any  unitary  rotation  Q.  leading  to  an  equivalent  stochastic  representation  for 
the  estimated  data  covariance  matrix: 


X  =  Ri/2q"z  RT1  =  (XX")  1  =  R-^Q"  (zz")  '  QR  l/2. 
Let  V  =  [v(0a)|v(0/f)],  and  choose  the  unitary  rotation  Q  such  that 


QR~1/2V  = 


^2x2 

°(/V-2)x2 


and  let  A  =  [<Sa|<5ft]. 


(A-7) 


(A-8) 


Note  from  (A-8)  that 


A"  A  =  V"R  *V  = 


Partition  the  white  data  matrix  such  that 

Z  = 


'  vH (0a)R~1\(0a) 

v"(0a)R'v(0,,)  ‘ 

’  ll^all2 

tab  ' 

_  vH(6>,,)R-1v(0„) 

vH(Gh)R  lv(9h)  _ 

1 - 

0 

IN|2  _ 

[^]{N-2)xL 


(A-9) 


(A-10) 


Using  well  known  matrix  inversion  lemmas  for  partitioned  matrices,  note  that  the  inverse  of  the 
whitened  sample  covariance  matrix  can  be  written  as 


(: ZZ")-1  = 


©2*2  * 


where  ©  =  Z 


1L  -  z2  (z^z2)  1  z^ 


z,. 


(a-ii: 


Note  that  0  is  white  complex  Wishart  distributed  with  L  —  N  +  2  degrees  of  freedom  [34].  denoted 
by  0  ~  CW(I2.  L  —  N  +  2).  Combining  this  matrix  partition  with  (A-8),  it  follows  that 

vH{0k)R~lv(ek)  =  vH(0k)R-1/2QH(ZZH)-lQR-1'2v(dk)  =  6%&-ldk  (A-12) 

for  k  =  a,  b.  Thus,  the  first  equation  in  (A-3)  can  be  written  as 


Pr  {Fa<F)  =  Pr  (SfJ &  lSb  —  F  ■  S  lSa  <  0) 

=  Pr  {tr  [0-1  (6bi 5'1  -  F  ■  6a6 ")]  <  ()} 
=  Prjtrf©-1  (Q,  A(F)Q")]  <()} 

=  Pr  {tr  [©_1  A(F)]  <  0}  , 


(A- 13) 
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The  last  equality  in  (A-13)  follows  from  the  invariance  of  complex  white  Gaussian  data  to  unitary 
rotations  [34,  49].  Define  the  following  variables  for  0: 


0  = 


’  01,1 

01,2  ' 

^0"’  = 

'  01’1 

0 1-2  ‘ 

1 

* — 

to 

to 

-01,2  ' 

02,1 

02,2 

02,1 

02,2 

=  ]e|  ’ 

-02,1 

01,1 

(A-14) 


From  (A-13)  it  follows  that 


fPr 

[0U  ^  — a2(f)1 

6>2’2  -  Ai(F) 

lPr 

01,1  .  -a2 (f)= 

02-2  '  A,(F) 

,  sign  [A  1(F)]  =  1 
,  sign[Ai(F)]  =  -1. 


From  matrix  inversion  lemma  in  (A-14)  note  that 


0  '  _  0^2  £  lXL-N+2 

Q2'2  01,1  2Xl-N+2 


(A-15) 


(A-16) 


where  the  last  equality  in  distribution  follows  from  the  fact  that  0  ~  CW(I->,L  —  N  +  2).  Thus, 
letting  l\(F)  =  — A2(F)/Ai(F),  it  follows  from  the  origin  of  (A-5)  and  (A-15)  that 

ri\(F) 

Pr  (FA  <  F)  =  0.5  •  {1  -  sign  [Ax(F)]}  +  sign[Ax(F)]  /  PFL-N+2tL^N+2(a)da.  (A-17) 

— OO 

Differentiating  (A-17),  the  desired  pdf  of  Fa  is  obtained: 


Pl\(F;  A)  = 


sign[At(F)]  •  PFl.n+2,l.n+2  [h(F)}  ■  —lx (F),  F  >  0 


0, 


otherwise 


(A-18) 


where 


i_  <L  \  ( ™  >  A^(F)  jL\  (tr\ 

dF  x  *  Ai (F)'dFX^  ^  +  [Aj(F)]2  '  dFAl^  ^ 


(A-19) 


Exact  closed  form  expressions  for  the  eigenvalues  Xn(F),  n  =  1,2  and  their  derivatives  are  derived  in 
Appendix  B.  Expressing  these  eigenvalues  in  terms  of  the  original  system  parameters  v(0a),  v(0^), 
and  R  is  accomplished  with  (A-9)  yielding  their  ratio  as  expressed  in  (29).  The  cdf  can  be  written 
compactly  as 


Pr  (Fa  <  F)  =  0.5  •  {1  —  sign  [Ai(F)]}  +  sign[Ai(F)]  •  F[l\(F),  L  -  N  +  2] 
from  which  the  pairwise  error  probabilities  easily  follow.  ■ 


(A-20) 
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A. 2  Bartlett  Algorithm 


Note  that  the  data  matrix  X  can  be  represented  as  a  linearly  transformed  white  data  matrix 
Z  that  is  invariant  to  any  unitary  rotation  Q,  leading  to  an  equivalent  stochastic  representation  for 
the  estimated  data  covariance  matrix: 

X  =  Ri/2q"z=>R  =  XXh  =  Kl/2qHZZHQR1/2.  (A-21) 

Let  V  =  [v(0a)|v(0/,)]  and  partition  the  whitened  data  matrix  as  in  (A-10).  Choose  the  unitary 
rotation  Q  such  that 


qrV^v  = 


A‘2x2 

°(Af-2)x2 


and  let  A  =  [5a|<5b]. 


Note  from  (A-22)  that 


A"A  =  v^rv  = 


'  v'V„)Rv(^) 

vH(Ga)Rv(Gb)  ' 

’  iNlI* 

6%  ' 

v"  (G/^Rvid,,) 

vH(dh) Rv{eh)  _ 

S^Sa 

INI2  . 

Note  that  the  whitened  sample  covariance  matrix  can  be  written  as 


ZZH  = 


(Z"Z,)2x2  * 


(A-22) 


(A-23) 


(A-24) 


Clearly,  the  matrix  Z{7  Z i  is  complex  central  Wishart  distributed  according  to  CW  >  (I2,  L)  [34,  49]. 
Combining  this  matrix  partition  with  (A-22),  it  follows  that 


vH(Gk)Rv(0k)  =  vH(ek)R^2Q"zz"QRl/2x{ek)  =  s'f!z[1zisk 

for  k  =  a,  b.  Thus,  the  second  equation  in  (A-3)  can  be  written  as 

Pr  (fa<f)  =  Pr  (d" Z" Z\S„  -  F  ■  Z[> Z , Sb  <  o) 

=  Pr  jtr  [Zf  Z,  (Sj"  -  F  •  Mfc")]  <  o} 
ZfZ,  (Q#A(F)Qg)]  <()} 


=  Pr  j  tr 
=  Pr  { tr 


Z"Z,A(F)]  <()}. 


(A-25) 


(A-2(>) 


The  last  equality  in  (A-26)  follows  from  the  invariance  of  complex  white  Gaussian  data  to  unitary 
rotations  [34,  49].  It  follows  from  (A-26)  that 


Pr 


ixl-MF)  +  2xl-HF)<0 


Pr 


sign  [A  1  (F)]  =  1 

ixj  < 

2Xl  _  > 

sign[Ai(F)]  =  -1 


-A  2(F) 
Ai(F) 


(A-27) 


Thus,  letting  lx(F)  =  —\‘2(F)/\i(F),  it  follows  from  the  origin  of  (A-5)  and  (A-27)  that 


ri-  ( F) 

Pr  (FA<F)  =  0.5  •  |l  —  sign  [a^F)]  |  +sign[A1(F)]  J  PFL,L(a)da. 
Differentiating  (A-28),  the  desired  pdf  of  F ^  is  obtained: 

p.  (F;  A)  =  (  sisn[A,{/-)J  ■  Pfl  l  [fj(F)]  •  ii^F),  F  >  0 

otherwise 


(A-28) 


(A-29) 


where 


0, 


-i l-x(F)  =  - iA 2(F)  +  } - iAj(F). 

dF  A  Ai(F)  dF  [Aj (F)]2  dF 


(A-30) 


Exact  closed  form  expressions  for  the  eigenvalues  Xn(F ),  n  =  L  2  and  their  derivatives  are  derived  in 
Appendix  B.  Expressing  these  eigenvalues  in  terms  of  the  original  system  parameters  v(0a).  v(<9/,), 
and  R  is  accomplished  with  (A-23).  The  cdf  can  be  written  compactly  as 


Pr(FA<F)  =  0.5  •  jl  —  sign  ^Aj(F)]  }  +  sign[Aj(F)]  •  F[/A(F),  L] 


(A-31) 


from  which  the  pairwise  error  probabilities  easily  followr.  ■ 


APPENDIX  B 

EIGENVALUES  OF  RANK  TWO  MATRICES 


It  is  desired  to  derive  closed  form  expressions  for  the  eigenvalues  and  their  derivative  with 
respect  to  F  of  the  following  parameterized  rank  two  matrix: 

A  =  aaH  —  F  •  bb,y  (B-l) 

where  clearly  a  and  b  are  linearly  independent.  Note  that  both  eigenvectors  must  be  linear  com¬ 
binations  of  these  two  vectors  such  that  if  Av  =  Av,  then  there  exists  some  complex  scalar  (3  such 
that  v  =  (3  •  a  +  b.  Satisfying  both  conditions  requires  that 


Av  =  Av 

(aa/y  —  F  •  bb/y)  •  ([3a.  +  b)  =  A  •  (/3a  +  b). 


(B-2) 


This  leads  to  two  relations  for  A 

A  =  ||a||2  +  (l//?)a"b 

A  =  -  (F/^a  +  FUbll2) 

that  when  satisfied  simultaneously  lead  to  the  following  quadratic  equation  for  the  scalar  (3 

(32  ■  (Fb/ya)  +  p  ■  (||a||2  +  F||b||2)  +  (a"b)  =  0.  (B-4) 


The  quadratic  formula  leads  to  the  following  solutions  for  the  eigenvalues 


Al.2(F) 


[Hall2-  ONI2! 

±J 

l|a||2  +  F||b||2' 

2 

2 

V 

2 

-  F\a"b\2 


that  can  be  readily  differentiated  to  obtain 


(B-5) 


(B-(i) 


APPENDIX  C 

COMPLEX  GAUSSI AN-BASED  STATISTICS:  A  SUMMARY 


This  appendix  provides  a  short  summary  of  statistics  related  to  the  complex  Gaussian  dis¬ 
tribution.  Similar  results  for  the  real  Gaussian  are  well  known  [34].  These  complex  analogs  are 
derived  in  [27]. 

Let  two  random  vectors  be  complex  normal  distributed  asz~  CM yv(0. 1/v)  and  is  ~  CAf/v(m,  I/v). 


Proposition  C.l  The  pdf  of  the  random  variable  p  =  ||z||2  is  given  by  a  complex  central  chi-square 
of  N  degrees  of  freedom: 

pN~le~'p/(N  -  1)!  p>  0.  (C-l) 

This  distribution  is  denoted  by  p  ~  \2N. 


Proposition  C.2  The  pdf  of  the  random  variable  ps  =  ||z<5||2  is  given  by  a  complex  noncentral 
chi-square  of  N  degrees  of  freedom: 

tf-'e-Ke-^oF^N^2  p5)/(N  -  1)!,  ps  >  0  (C-2) 

where  oFi(-)  is  a  hypergeometric  function  [ 1 /,  and  the  noncentrality  parameter  is  given  by  S2  = 

1 1 rri 1 1 2 .  This  distribution  is  denoted  by  ps  ~  XjvW* 


Proposition  C.3  The  pdf  of  the  random  variable  F^,a/  =  where  both  chi-squared  statistics 

are  independent  is  given  by 


(N  +  A/-1)!  F"-1  ~> 

(N  -  l)!(iV/  -  1)!  (l  +  F)(w+A/)’  ~  ’ 


(C-3) 


a  complex  central  F -distribution  with  N,M  degrees  of  freedom. 


Proposition  C.4  The  pdf  of  the  random  variable  F/v,A/(^)  =  XatW/Xa/  where  both  chi-squared 
statistics  are  independent  is  given  by 


{N  +  M-l)\  fn_1 
(N  -  1)!(M  -  1)!  (l  +  F)(jv+a/) 


e~s%Fi[N  +  M;N;62  F/(l  +  F)],  F  >  0 


(C-4) 


a  complex  noncentral  F- distribution  with  N ,  M  degrees  of  freedom,  where  \F\( •;  •;  •)  is  the  confluent 
hypergeometric  function  [1].  The  cumulative  distribution  function  (cdf)  for  this  random  variable  is 
given  by 


Fr(FN,A/  <  x)  = 


(1  +x) 


N+M 


A/-1 

-E 

k=0 


N  +  M  -  1 
h  +  N 


IGk 


+  i 


X2 


1  T  x 


(C-5) 


where  IGm{a)  =  e  a  ak/h!  is  the  incomplete  Gamma  function.  The  cdf  for  the  central 

complex  F  is  the  special  case  of  6  =  0. 
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Proposition  C.5  The  pdf  of  the  random  variable  Pn.m  =  1/(1  +  Fa/,jv)  is  given  by 


(N  +  M- 1)!  ~N_, 
(N-  1)!(M-  1)!' 


0  <  0  <  1 


(C-6) 


a  complex  central  beta  distribution.  In  addition,  its  cdf  is  given  by  the  convenient  finite  t 


sum 


(C-7) 


Noting  that  Pr(F/v,A/  <  x)  =  Pr(— 1  +  l /Pa /,yv  <  x),  a  second  representation  for  the  cdf  of  the 
central  F  is  possible. 


Proposition  C.6  The  pdf  of  the  random  variable  $n,m(8)  =  1/(1  +  Fm,n(8))  is  given  by 

_{N  +  M- 1)1  ~N_i  _  ~  M_x  '  y _ N\(M  —  !)! _ (^T  _  fan 

(N  —  1)!(A/  —  1)!  1  P>  ^  (N  -  n)\(M  -  1  -  »)!  n!  1  P)  ' 


72=0 


(C-8) 


for  0  <  P  <  1.  is  known  as  the  complex  noncentral  beta  distribution. 

Proposition  C.7  The  cdf  of  the  random  variable  Pn,m(8)  =  1/(1  +  Pa/,/v(^))  is  given  by 

N  4-  A/ — 1  /  _  r  ,  ,  l\/l  \  777, 

<  x)  =  1  -  x^-1  £  I  +  J  ~  (“ T  /  IGm-M+i  {x  ■  62)  (C-9) 

772  =  A/  \  /  '  *  ' 

/or  a;  >  0. 
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APPENDIX  D 

MATLAB  CODE  FOR  CANONICAL  EXAMPLE:  SINGLE  BROADSIDE 
PLANEWAVE  SIGNAL  IN  WHITE  NOISE 


1  function  []  =  Capon. Alg.MSE; 

2  •/. - 

3  7.  This  code  computes  the  threshold  region  mean  squared  error  (MSE)  performance  of 

4  '/,  the  Direction-Of -Arrival  (DOA)  estimate  obtained  from  the  Capon-MVDR  algorithm  of 

5  7.  a  single  planewave  signal  in  white  noise  as  observed  on  a  uniform  liinear  array  (ULA)  . 

6  '/,  This  code  is  based  upon  the  theoretical  results  to  be  published  in  the  IEEE 

7  '/.  Transactions  on  Signal  Processing,  2005. 

8  7. 

9  7,  Christ  D.  Richmond,  christQll.mit.edu,  12/10/2004 

10  7,  MIT  Lincoln  Laboratory,  All  rights  reserved 

11  7. - 

12  7.  SET  RUN  PARAMETERS.  .  . 

13  7. - 

14  7.  Define  GLOBAL  variables 

15  global  run.params; 

16 

17  7.  Define  model  parameters .  .  . 

18 

19  run.params. N  =  18;  7t  Number  of  array  sensors,  beams,  or  spatial  channels 

20  run.params. d  =  1/2.25;  7.  Element  separation  in  wavelengths 

21  run.params. L  =  3  *  run.params .  N ;  7.  Sample  support 

22  run.params . K  =  run.params .  L  -  run.params. N  +  2;  7*  Degrees  of  freedom  for  F-distribution 

23  7.  of  pairwise  error  probability 

24  run.params .  wl_dB  =  0;  7.  White  noise  level  dB 

25  run.params .  Ang.low  =  90-30;  7.  Lower  limit  of  search  angles  (degrees); 

26  run_params.Ang.high  =  90+30;  7.  Upper  limit  of  search  angles  (degrees); 

27  run.params  .Del.Ang  =  run.params  .  Ang.high  -  run.params  .  Ang.low;  7.  Degrees 

28  7. - 

29  7.  MVDR  /  Capon  Parameters .  .  . 

30  7. - 

31  run.params . Sig.angles  =  [90];  7.  True  signal  direction  (degrees;  90degs  is  broadside) 

32  see_src.no  =  1;  7.  Source  #  whose  MSE  is  of  interest,  for  multiple  signal  case 

33  7.  [not  included  here] 

34  run.params . Tgt.angs.degs  =  run.params . Sig.angles ; 

35  run_p2Lrams.No.sigs  =  length  (run.params  .  Sig.angles)  ; 

36  7. - 

37  run.params  .  SNRs.dB  =  1  inspace (-30, 15, 50) -run.params  .wl.dB;  7.  Element  level 

38  7.  signal-to-noise 

39  7.  ratio  (SNR)  (decibels) 

40  run.params  .  SNRs  =  bd (run.params  .  SNRs.dB)  ;  7.  Element  level  SNR  in  power  domain 

41  run_paLrams.ULA_BW.degs  =  ULA_3dB_BW(run_p2Lrams  .N, run.params  .d)  ;  7.  Beamwidth  for  ULA  (degs) 

42  7. - 

43  7.  GET  TRUE  DATA  COVARIANCES .  .  . 

44  7. - 

45  R.mvdr  =  cell(l , length (run.params . SNRs) ) ; 

46  Ri.mvdr  =  cell (1 , length (run.params . SNRs) ) ; 

47  R.mvdr _sqrt  =  cell (1 , length (run.params . SNRs) ) ; 

48  R.mvdr.si  =  cell (1 , length (run.params . SNRs) ) ; 
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49 

50 

51 

52 

53 

54 

55 

56 

57 

58 

59 
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62 

63 

64 

65 

66 

67 

68 

69 

70 

71 

72 

73 

74 

75 

76 

77 

78 

79 

80 
81 
82 

83 

84 

85 

86 

87 

88 

89 

90 

91 

92 

93 

94 

95 

96 

97 

98 

99 

100 
101 
102 


for  usnr  =  1 : length (run_params . SNRs) ; 

'/,  function [R,Ri]  =getRvan(N,wl_dB,tht ,thl_dB,d,f lg) ; 

[R,Ri]  =  getRvan(run_params . N ,run_params . wl_dB , [run.params . Sig.angles] ,  ... 

[run.params . SNRs.dB (usnr) *ones (1 ,run_params . No.sigs)]  , run.params . d) ; 
R_mvdr{usnr}  =  R;  Ri_mvdr{usnr}  =  Ri; 

'/,  function  [Rsqrt  ,Rsqrti  ,DetRsqrt]  =  get.Rsqrt (R) ; 

[Rsqrt ,Rsqrti, jnk]  =  get.Rsqrt (R) ;  R_mvdr_sqrt{usnr}  =  Rsqrt;  R_mvdr_si{usnr}  =  Rsqrti 

end;  '/,  for  usnr  =  1 :  length (run_params  .  SNRs)  ; 

run_params . R_mvdr  =  R_mvdr;  run_params . Ri_mvdr  =  Ri_mvdr; 

run.params . R_mvdr_sqrt  =  R_mvdr_sqrt;  run_params . R_mvdr_si  =  R_mvdr_si; 

•/. - 

run_p2LT2Lms  .res_f ac  =  40;  '/,  Conventional  beamsplit  ratio  for  angle  search 

'/.  #  grid  samples  in  angle  search.  .  . 

run_params . Nres  =  run_params . res_f ac*ceil(180/run_params . ULA.BW.degs) ; 
anggs  =  linspace (run.params . Ang.low , run.paxams . Ang_high , run.params . Nres) ; 

V  =  van(run_params  .N, anggs  ,run_params  .d)  ;  '/,  search  angles  and  steering  vectors 
run.params  .  Vanggs  =  van  (run_params  .  N  ,  anggs  ,run_params  .  d)  ;  run_params  .  anggs  =  anggs; 

*/. - 

run_params  .  angsi  =  run_parauns  .  Sig_angles  ; 

run.paLrams  .  Vs  =  van  (run _par  am s  .  N ,  run_params  .  angsi  ,  run_params  .  d) ; 
run_pairaLms  .  V  =  van  (run.parauns  .  N ,  run.params  .  anggs  ,  run_params  .  d)  ; 

7, - 

y.  COMPUTE  ASYMPTOTIC  LOCAL  ERROR  MSE  AND  BOUNDS .  .  . 

7, - 

'/,  Get  Cramer-Rao  Bound.  .  . 


*/,  See  C.  D.  Richmond,  "Craimer-Rao  Bounds:  The  Adaptive  Array  Story,"  Proceedings  of 
*/,  the  Adaptive  Sensor  Array  Processing  Workshop,  MIT  Lincoln  Laboratory,  March  2000. 

disp(  [  Computing  CRB  Local  Bounds...’]); 


thT  =  run.params . Tgt_angs_degs ; 

VT  =  vanCrun.paxams  .  N , thT,run_paxams  .  d)  ;  '/»  True  target  array  response 
den_ang  =  zeros (length (run_paxams . SNRs) , length (run.params .  angsi) ) ; 

for  usnr=l : length (run.params . SNRs) ;  u  =  usnr;  sig_S  =  run.params . SNRs (usnr) ; 
run.params.R  =  run.paLrams . R_mvdr{usnr} ;  run.params . Ri  =  run_params . Ri_mvdr{usnr} ; 

x - 

D  =  van_p(run_params.N,thT,run_params .d) ; 

dRsdw  =  hm(  sig_S  *  (  D*VT>  +  VT*D’  )  )  ;  dRsdwds  =  hm(  D*VT’  +  VT*D’  ); 
dRsds  =  hm(  VT^VT’  ); 

Rs  =  run_params  .  R_mvdr{usnr> ;  Rsi  =  run_paurams  .  Ri_mvdr{usnr} ; 


F  =  zeros(2);  F(l,l)  =  trace( 
F(l,2)  =  trace (  dRsdw  *  Rsi  * 
F(2,2)  =  trace (  dRsds  *  Rsi  * 
•/. - 


dRsdw  *  Rsi  *  dRsdw  *  Rsi  ) ; 
dRsds  *  Rsi  ) ; 

dRsds  *  Rsi  );  F(2 , l)=conj (F(l , 2) ) ; 


F  =  run_params . L  *hm(F);  J=hm(  inv(F)  ); 
crbw(usnr)  =  l/abs(  F(l,l)  );  crbws(usnr)  =  abs(  J(l,l)  ); 
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103  den_ang(usnr , : )  =  l/abs(  J ( 1 , 1 )  ); 

104  end;  */,  for  usnr=l : length (run.params . SNRs) ;  u  =  usnr; 

105 

106  crb  =  l./den_ang;  '/,  CRB  in  radians 

107  X - 

108  '/,  Get  Vaidyanathan-Buckley  (VB)  local  error  Capon  MSE  .  .  . 

109  '/,  Based  on  C.  Vaidyanathan,  K.  M.  Buckley,  “Performance  Analysis  of  the  MVDR  Spatial 

110  '/,  Spectral  Estimator,”  IEEE  Transactions  on  Signal  Processing,  Vol.  43,  No.  6, 

111  7.  pp.  1427—1437,  June  1995. 

112 

113  thT  =  run_params . Tgt_angs_degs ; 

114  VT  =  van(run_params .  N,thT,run_params  .d)  ;  */,  True  target  array  response 

115  den_ang  =  zeros ( length (run_params . SNRs) , length (run_params . angsi) ) ; 

116 

117  for  usnr=l : length (run_params . SNRs) ; 

118 

119  run.params.R  =  run_params . R_mvdr{usnr} ;  run_params . Ri  =  run.params . Ri_mvdr{usnr} ; 

120  Rs  =  run.params . R_mvdr{usnr> ;  Rsi  =  run.params . Ri_mvdr{usnr> ; 

121  thT.true  =  run_params . Sig.angles (see_src_no) ; 

122  */. - 

123  '/,  Get  asymptotic  bias... 

124  a  =  van(run_params . N , thT.true , run.params .d) ; 

125  ap  =  van_p(run_params . N , thT.true ,run_params . d) ; 

126  app  =  van_p2(run_params.N,thT_true,run_params.d) ; 

127  Ri  =  run_params .Ri_mvdr{usnr}; 

128 

129  thT.asymp  =  -  real(  ap’*Ri*a  )  /  real(  ap’*Ri*ap  +  a’*Ri*app  )  +  thT.true; 

130  Asymp_Bias  =  thT_asymp  -  thT_true;  '/,  Only  accurate  when  bias  is  small... 

131  7. - 

132  */.  Evaluate  at  ASYMPTOTIC  MSE  of  angle  estimates... 

133  a  =  van(run_params .N,thT_asymp,run_params .d) ; 

134  ap  =  van_p(run_params . N ,thT_asymp,run_params . d) ; 

135  app  =  van_p2(run_params . N,thT_asymp,run_params .d) ; 

136  ap3  =  van_p3(run_params  .  N,thT_asymp,run_pairams  .d)  ; 

137 

138  '/,  use  notation  of  VB  paper  for  convenience... 

139  L  =  run.params . L ;  N  =  run_params . L ;  K  =  run.params . N ; 

140 

141  fpp  =  2  *  real(  ap’*Ri*ap  +  a’*Ri*app  );  B  =  ap*a’  +  a*ap’; 

142  fp3  =  2  *  real(  3  *  ap’*Ri*app  +  a’*Ri*ap3  ); 

143  ED 2  =  (2  *  (N-  K)  /  (  (N-K-l) * (N-K+l)  )  )  *  real(  a’*Ri  *B*Ri*ap)/(  fpp"2  ); 

144  ED3  =  (2  *  (N-  K)  /  (  (N-K-l) * (N-K+l)  )  )  *  real(  a’*Ri  *  B  *  Ri  *  app  )  /  (  fpp~2  ); 

145  AddBias  =  ED3  -  fp3  *  ED2  /  (  2  *  fpp  ); 

146 

147  '/,  Evaluate  using  BOTH 

148  Bias  =  Asymp.Bias  +  AddBias;  den_ang(usnr ,  : )  =  l/(Bias~2  +  ED2  );  */,  1/(VB  MSE)  in  1/radians 

149  end;  '/,  for  usnr=l :  length (run.params  .  SNRs) ; 

150  % - 

151  y.  COMPUTE  GLOBAL  ERRORS  USING  METHOD  OF  INTERVAL  ERRORS  (MIE)  .  .  . 

152  X - 

153  disp( [ ’Determining  location  of  local  maxima/peaks/sidelobes ...’]) ; 

154 

155  peak_SLL_idx  =  cell (length (run.params . angsi) , 1) ; 

156  peak_SLL_angs  =  cell (length (run_params . angsi) , 1) ; 
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157  peak.SLL.V  =  cell (length (run.params . angsi) , 1) ; 

158  peak_SLL_vals  =  cell (length (run.params . angsi) , 1) ; 

159 

160  ridx  =  length ( run.params . SNRs.dB) ; 

161  y. - 

162  mx_dim  =  0;  usnr  =  length (run.params . SNRs) ;  ua  =  see_src_no; 

163  Vo  =  run_params . V ;  Ri  =  run_params . Ri_mvdr{usnr} ; 

164  t.mvdr  =  l./sum(  conj (  Vo  )  . *  (  Ri*Vo  )  ,  1  );  tMVDR.dB  =  dbl0(  t.mvdr  ); 

165 

166  f.ambigO  =  figure;  plot(  run.params . anggs ,  tMVDR.dB  ); 

167  xlabel( ’Search  Angle  (degrees)’);  ylabel(’SNR  (dB)’); 

168  title (’MVDR  Ambiguity  Function’);  vln (run.params . Sig.angles , ’r- .’) ; 

169 

170  irm.idx  =  f ind(imregionalmax(tMVDR_dB) ) ; 

171 

172  '/,  Throw  out  maxima  due  to  true  signals.  .  . 

173  idx  =  [] ;  for  u  =  1 : run.params . No_sigs ; 

174  idx  =  [idx  find(abs (anggs (irm.idx) -run.params. Sig.angles (u) ) <=run_params .ULA_BW_degs)] ; 

175  end;  '/,  for  u  =  1 : run.params  . No.sigs ; 

176 

177  irm.idx  =  irm_idx(  setdif f ( [1 : length (irm.idx)] , idx)  ); 

178  aa  =  tMVDR.dB (irm.idx) ;  [jnk.iro]  =  sort(  aa  );  iro  =  flipud(  iro(:)  ).’; 

179  f igure (f .ambigO) ; 

180  hold  on;  plot (run.params . anggs (irm.idx) , tMVDR.dB (irm.idx) , ’ko ’) ;  hold  off; 

181 

182  '/,  If  local  max  is  not  an  ENDPOINT,  then  refine.  .  . 

183  run.params .Ri.hat  =  run.params . Ri_mvdr{usnr} ; 

184  for  u.nsig  =  1 : length(irm_idx(iro) ) ; 

185  '/.  refine  if  NOT  endpoint  max.  .  . 

186  if  (  irm.idx (iro (u.nsig) )  ~  =  1  )  &  (  irm.idx (iro (u.nsig) )  ~=  length (tMVDR.dB)  ); 

187  idx  =  irm.idx (iro (u.nsig) ) ; 

188  fh  =  QMVDR.fnc;  */,  Use  function  handle  to  make  visible  to  fminbnd  (a  MATLAB  function); 

189  a.refine (u.nsig)  =  ... 

190  f minbnd(fh, max (run.params . Ang.low, anggs (idx) -(1/2) *run_params . ULA.BW.degs) , . . . 

191  min (anggs (idx) +(1/2) ♦run.params . ULA_BW_degs .run.params . Ang.high) ) ; 

192 

193  '/,  Check  that  angle  estimate  refinement  is  small... if  not,  then  don’t  use... 

194  7,  could  be  a  saddle  point.  .  . 

195  if  abs(  anggs (idx)  -  a.refine (u.nsig)  )  >=  ( 1/4) ^run.params .ULA.BW.degs ; 

196  a.refine (u.nsig)  =  anggs(idx);  end; 

197 

198  else; 

199  a.refine (u.nsig)  =  anggs(  irm.idx (iro (u.nsig) )  ); 

200  end;  '/,  if  (  irm.idx  (iro  (u.nsig) )  ~=  1  )  &  (  irm.idx  (iro  (u.nsig) )  ~=  length  (tMVDR.dB)  ); 

201 

202  peak.SLL.angsfua} (u.nsig)  =  a.ref ine (u.nsig) ; 

203  peak_SLL_V{ua}( : , u.nsig)  =  van( run.params . N ,peak_SLL_angs{ua} (u.nsig) , run.params . d) ; 

204 

205  end;  '/,  for  u.nsig  =  1 :  length  (irm.idx)  ; 

206 

207  aV  =  peak_SLL_V{ua} ;  peak_SLL_vals{ua}  =  dbl0(  l./abs(  sum(  conj (  aV  )  .*  (Ri*aV) ,  1)  )  ); 

208  aa  =  dbl0(  l./abs(  sum(  con j (  aV  )  .*  (Ri*aV) ,  1)  )  ); 

209  if  length(aa)  >  mx.dim;  mx.dim  =  length ( aa) ;  end; 

210  run.params . peak.SLL.angs  =  peak.SLL.angs ;  run.params . peak.SLL.V  =  peak.SLL.V; 


211  peak_SLL_idx  =  irm_idx(iro) ; 

212  y, - 

213  7.  BUILD  MSE  ERROR  TERMS  FOR  GLOBAL  ERROR  CONTRIBUTIONS... 

214  7. - 

215  disp( ’Building  Error  and  mu/lambda  Matrices...’); 

216  Err.Mat  =  zeros (1 ,mx_dim) ;  Ptht_cond  =  zeros (1 ,mx_dim) ;  ua  =  see_src_no; 

217  possible_angs  =  peak_SLL_angs{ua> ; 

218  for  upa  =  1 : length (possible_angs) ; 

219  Err_Mat (1 ,upa)  =  (pi/180)~2  *  abs(  run.params . angsi (ua)  -  possible_angs (upa)  )~2;  end; 

220  7. - 

221  7.  COMPUTE  TOTAL  MSE  USING  MIE  MSE  APPROXIMATION.  .  . 

222  7. - 

223  MF_MSE_Approx  =  zeros (length (run.params . SNRs) , 1) ; 

224  Glb.MSE  =  zeros (length (run.params . SNRs) , 1) ;  Loc_MSE=  zeros (length (run.params . SNRs) , 1) ; 

225  tO  =  clock;  rand (’ state ’, sum (100*clock) ) ; 

226  disp( [ ’Beginning  MSE  Calculations...’]); 

227  7,  True  array  response 

228  run.params . vT  =  van(run_params . N ,run_params . Sig.angles (see_src_no) ,run_params . d) ; 

229 

230  for  usnr  =  1 :  length  (run.params .  SNRs)  ;  7*  Compute  MSE  as  function  of  element  level  SNR... 

231  Prob_Err_Mat  =  zeros (1 ,mx_dim) ; 

232  for  upa  =  1 :  length (possible_angs)  ;  7*  Errors 

233  7.  function  [Pe]  =  Pe_Capon(R, va, vb,L)  ; 

234  Prob_Err_Mat (1 ,upa)  =  ... 

235  Pe.Capon (run.params . Ri_mvdr {usnr} , peak_SLL_V{ua} ( : , upa) , run_params . vT , run.params . L) ; 

236  end;  7.  for  upa  =  1 :  length(possible_angs)  ; 

237 

238  MF_MSE_Approx(usnr)  =  ... 

239  (1 . /den.ang(usnr) ) . ’  .*  (l-min(l ,sum(Prob_Err_Mat ,2) ) )+sum(  Err.Mat  .*  Prob_Err_Mat ,  2  ); 

240  amidon ( tO , usnr , length (run_par ams . SNRs ) ) ; 

241  end;  7.  for  usnr  =  1 : length (run.params . SNRs) ; 

242 

243  run.params .MF_MSE_Approx  =  MF_MSE_Approx ;  MF_MSE_ApproxO  =  MF_MSE_ Approx ; 

244  7. - 

245  7.  Clip  MSE  so  that  it  maxes  out  at  MSE  of  UNIFORMLY  distirbuted  estimate... 

246 

247  BW.degs  =  run.params .ULA_BW_degs ; 

248 

249  xx  =  linspace (run.params . Ang.low, run.params . Ang.high, 1000) ; 

250  fxx  =  (  l/(  xx(end)  -  xx(l)  )  )  *  ones (size (xx) ) ; 

251  varx  =  trapz(  xx,  fxx  . *  (xx  -  run.params . Sig.angles (see_src.no) ). "2  ); 

252  unif.mse  =  varx  *  (pi/180) "2; 

253 

254  for  u=l : length (MF.MSE.Approx) ; 

255  MF.MSE.Approx  (u)  =  min  (MF.MSE.Approx  (u)  ,  unif.mse)  ;  7.  Clipping... 

256  end;  7.  for  u=l :  length  (MF.MSE.Approx)  ; 

257  7. - 

258  7.  PLOT  FINAL  RESULTS .  .  . 

259  f igure ; 

260  plot (run.params . SNRs.dB ,dbl0(  sqrt(  (180/pi) "2  *(  MF.MSE.Approx  )*... 

261  ( 1 /run.params . ULA.BW.degs) "2  )  ),’r-’); 

262  hold  on; 

263  plot (run.params . SNRs.dB ,dbl0(  sqrt(  (180/pi)~2  *(  crb  )*... 

264  ( 1 /run.params . ULA.BW.degs) ~2  )  ),’k:’); 


265  hold  off ; 

266  hold  on; 

267  plot (run.params . SNRs.dB ,dblO(  sqrt(  (180/pi)~2  *(  l./den_ang  )*... 

268  (l/run_params . ULA_BW_degs) "2  )  ),’g — ’); 

269  hold  off; 

270  titled  Capon  Algorithm  Threshold  Region  MSE  Prediction’); 

271  xlabel( ’Element  Level  SNR  (dB)’);  ylabel(’Root  MSE  in  Beamwidths  (dB)’) 

272 

273  '/,  Direct  threshold  SNR  approximation  for  Capon  algorithm.  .  . 

274  [jnk,ioi]  =  max(  abs(  bd(  peak_SLL_vals{l}  )  )  ); 

275  vl  =  run.params . vT ;  v2  =  peak_SLL_V{l}( : , ioi) ;  al2  =  abs(  vl’*v2  )~2; 

276  al  =  abs(vl’*vl);  a2  =  abs(v2’*v2);  co  =  al2/al/a2;  gam  =  al ; 

277  [Est_Thg]  =  Capon_ThSNR_II_dB(run_params . N, co ,run_params . L,run_params . N) ; 

278  vln(dblO(Est_Thg),’k-.’);  axis([-30  10  -20  5]); 

279  legend ( ’ MIE  MSE ’ , ’ CRB ’ , ’ VB  MSE ’ , ’ Approx .  SNR.TH ’ ) 

280  return; 

281  7. - 

282  7,  ALL  SUBFUNCTIONS .... 

283  7. - 

284  function  [Pe]  =  Pe_Capon(Ri , va, vb,L) ; 

285 

286  7.  function  [Pe]  =  Pe_Capon(Ri , va, vb,L)  ; 

287  7. 

288  7.  This  code  implements  the  algorithm  for  computing  the  pairwise  error  probability  for 

289  7.  the  Capon  algorithm  presented  in  IEEE  T-SP  paper  by  C.  D.  Richmond,  to  appear  2005. 

290  7. 

291  7.  INPUTS: 

292  7. - 

293  7.  Ri  -  Inverse  of  true  data  covariance  matrix 

294  7.  va  -  Array  response  for  error/outlier  sidelobe  direction 

295  7.  (location  of  local  max  of  ambiguity  fnc) 

296  7.  vb  -  Array  response  for  true  signal 

297  7.  L  -  Number  of  data  samples  (spatial  snapshots)  used  in  covariance  estimation 

298  7. 

299  7.  OUTPUT: 

300  7. - 

301  7.  Pe  -  Pairwise  error  probability  of  Capon  algorirth: 

302  7.  Pe  =  Prob[  P_Capon(tht_a)  >  P_Capon(tht_b)  ] 

303  7. 

304  7.  Author:  Christ  D.  Richmond 

305  7.  Date:  8/18/04 

306 

307  N  =  length (va) ;  K  =  L  -  N  +  2;  F  =  1 ; 

308  pa  =  abs(  va’*Ri*va  ) ;  pb  =  abs(  vb’*Ri*vb  );  lab  =  F  *  abs(  va’*Ri*vb  )~2; 

309  Pp  =  (pb  +  F  *  pa)/2;  Pm  =  (pb  -  F  *  pa)/2;  D  =  sqrt(  Pp~2  -  lab  ); 

310  1  =  -(  Pm  +  D  )  /  (Pm  -  D);  S  =  sign(Pm  -  D) ; 

311  if  K  <=  40; 

312  Pe  =  0.5  *  (1  +  S)  -  S  *  cumF(K,K,0,l) ; 

313  else; 

314  sigb  =  sqrt (  (K+l)/(2  *  (2*K  +  1  )  )  -  (0.5)~2  ); 

315  arg.ef  =  (  1/ ( 1+1)  -  0.5  )/(  sigb  *  sqrt(2)  ); 

316  Pe  =  0.5  *  (1  +  S)  -  S  *  0.5  *  erf c(arg_ef ) ; 

317  end;  7.  if  K  <=  40; 

318 
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319  return; 

320  % - 

321  function  [Th]  =  Capon_ThSNR_II_dB(gam, co ,L,N , eta) ; 

322 

323  7.  function  [Th]  =  Capon_ThSNR_dB(gam,co,L,N) ; 

324  7. 

325  7.  This  function  implements  the  direct  computation  of  an  approximation  for 

326  %  the  threshold  SNR  of  the  Capon  algorithm  for  the  fundamental  and  canonical  case  of 

327  7.  a  single  planewave  source  in  white  noise.  The  approximation  is  valid  in  the 

328  7.  large  sample  case,  i.e.  approximately  L  >=  2N. 

329  7. 

330  7.  INPUTS : 

331  7. - 

332  7.  gam  -  Norm  of  array  responses  (equal  norm  is  assumed)  .  Note  that  gam  =  #  sensors 

333  7.  for  usual  Vandermonde  steering  vectors. 

334  '/,  co  -  Geometric  cosine  between  true  signal  array  response  and  array  response  for 

335  7,  direction  of  highest  sidelobe  of  Capon  ambiguity  function. 

336  7.  L  -  Sample  support  of  estimated  data  covariance 

337  7.  N  -  number  of  sensors  in  array 

338  7. 

339  7.  OUTPUT: 

340  7. - 

341  7.  Th  -  Capon  threshold  SNR  (in  power  domain,  not  in  decibels) 

342  7. 

343  7.  Author:  Christ  D.  Richmond 

344  7.  Date:  8/18/2004 

345  7. - 

346  M  =  L  -  N  +  2;  n.std  =  4; 

347  if  nargin  <  5; 

348  eta  =  -  (  1  -  n.std  *  sqrt(  l/(  2*M  +1)))/(1+  n.std  *  sqrt(  l/(  2*M  +  1  )  )  ); 

349  end;  7.  if  nargin  <  5; 

350  7. - 

351  xi  =  (  (1-eta) /(1+eta)  )~2;  em  =  (xi  -  1);  ep  =  (xi  +  1);  psi  =  ep/em;  so  =  1-co; 

352  7. - 

353  7.  coefficients  of  quartic... 

354  a  =  2*gam*so  +  (2/gam) * (1-psi) ; 

355  b  =  (4*co/gam/gam/em)  +  2*(  (gam*gam  *  so*so  /  2)  -  co)  +  2*(l-psi)*(  1/gam/gam  +  1+so  ); 

356  c  =  (8*co/gam/em)  +  2*(  (l+so)/gam  +  gam*so*so  )  -  2*psi*(  so* (gam+l/gam)  +  1/gam); 

357  d  =  4*co/em  +  co~2  +  2*so* ( 1-psi) ; 

358  7. - 

359  7.  Apply  solution  to  quartic  equation.  .  . 

360  p  =  b;  q  =  (a*c  -  4*d) ;  r  =  (a~2  *  d  +  c~2  -  4*b*d) ;  u  =  q  -  p~2  /  3; 

361  v  =  r  -  p*q/3  +  2*p~3/27; 

362 

363  jO  =  4*(u/3)~3  +  v~2;  7.  Discriminant 

364 

365  if  jO  >  0; 

366 

367  w  =  sqrt  (  jO  )  ;  7.  y  =  (  (w  -  v)/2  )~(l/3)  -  (u/3)  *  (2/  (w-v) )  ~  (1/3)  -  p/3; 

368  y  =  (u/3) * (2/ (w+v) ) ~ (1/3)  -  ( (w+v) /2) - (1/3)  -  p/3; 

369 

370  elseif  jO  <=  0; 

371 

372  s  =  sqrt(  (-u/3)  );  t  =  -v/(2*s~3);  k  =  acos(t)/3; 


373  ck  =  cos(k);  sk  =  sin(k); 

374  yrts  =  real(  [  2  *s*ck  -  p/3;  s  *  (-ck  +  sqrt(3)  *  sk  )  -  p/3;  ... 

375  s  *  (-ck  -  sqrt(3)  *  sk  )  -  p/3]  ); 

376 

377  y  =  yrts(l);  7.  All  roots  appear  to  work... 

378 

379  end;  7,  if  jO  >  0; 

380 

381  e2  =  a~2/4  -  b  -  y;  e  =  sqrt(e2);  f2  =  (y . “2)* (1/4)  -  d; 

382  f  =  sqrt(f2);  ef  =  a*y*(l/4)  +  c/2; 

383  G  =  a/2  +  e;  g  =  a/2  -e;  H  =  -y./2  +  f;  h  =  -y./2  -  f; 

384  frts  =  [  (  -G  +  sqrt (  G~2  -  4*H  )  )/2;  (  -G  -  sqrt(  G~2  -  4*H  )  )/2; 

385  (  -g  +  sqrt (  g~2  -  4*h  )  )/2  ;  (  -g  -  sqrt(  g~2  -  4*h  )  )/2] ; 

386  Th  =  frts  (3);  7.  Root  leading  to  desired  solution 

387  return; 

388  7. - 

389  function [x] =bd(y) ;  x  =  10 . “ (y . /10) ;  return; 

390  7. - 

391  function [y] =dbl0(x) ;  y  =  10*logl0(  abs(x)  );  return; 

392  7. - 

393  function  [H]  =  hm(  A);H=(A+A’  )  *  (1/2);  return; 

394  7. - 

395  function  [out]  =  ULA_3dB_BW(N ,d) ; 

396 

397  7.  function  [out]  =  ULA_3dB_BW(N  ,d)  ; 

398  7. 

399  7.  Computes  beamwidth  for  uniform  linear  array. 

400  7. 

401  7.  N  -  number  of  elements 

402  7.  d  -  element  separation  in  wavelengths 

403  7. 

404  L  =  N*d;  out  =  (180/pi)  *  (1  ./  L);  return; 

405  7. - 

406  f  unct ion [v] =van (N , ang , elem , f lg) 

407 

408  7.  f  unct  ion  [v]=van(N,  ang,  elem) 

409  7. 

410  7.  This  function  computes  simple  Vandermonde  steering  vectors. 

411  7. 

412  7.  N  -  Number  of  elements 

413  7.  ang  -  Look  direction  in  degs :  0  deg  endfire;  90  deg  broadside 

414  7.  180  opposite  endfire 

415  7.  elem  -  Element  spacing  in  wavelenths 

416  7#  fig  -  Flag  to  make  array  center  the  phase  reference 

417 

418  if  nargin  <  4;  v  =  exp( j*2*pi*elem* [0 : N-l] ’ *cos (ang( : ) ’ *pi/180) ) ; 

419  else;  v  =  exp( j*2*pi*elem* (  [0 : N-l] *  -  (N-l)/2  ) *cos (ang( : ) ’ *pi/180) ) ;  end; 

420  7. - 

421  f  unct ion [v] =van_p (N , ang , elem , f lg) 

422 

423  7.  f  unct  ion  [v]  =van_p(N,  ang,  elem,  fig) 

424  7. 

425  7.  This  function  computes  the  first  derivative  of  Vandermonde  steering  vectors 

426  7. 
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427  7.  N  -  Number  of  elements 

428  7.  ang  -  Look  direction  in  degs:  0  deg  endfire;  90  deg  broadside 

429  7,  180  opposite  endfire 

430  7.  elem  -  Element  spacing  in  wavelenths 

431 

432  if  nargin  <  4; 

433  C  =  - j*2*pi*elem* [0:N-1] ’*sin(ang( : ) ’*pi/180) ; 

434  v  =  C  . *  exp(j*2*pi*elem* [0 : N— 1] ’ *cos(ang( : ) ’ *pi/180) ) ; 

435  else; 

436  C  =  - j*2*pi*elem* ( [0:N-1] ’  -  (N-l)/2  )*sin(ang( : ) 7 *pi/180) ; 

437  v  =  C  .*  exp(j*2*pi*elem*( [0 : N— 1] ’  -  (N-l)/2  ) *cos (ang( : ) ’ *pi/180) ) ; 

438  end; 

439  return ; 

440  7. - 

441  function [v] =van_p2(N, ang, elem, fig) 

442 

443  7.  function [v] =van_p2(N, ang, elem, fig) 

444  7. 

445  7.  This  function  computes  the  second  derivative  of  Vandermonde  steering  vectors. 

446  7. 

447  7.  N  -  Number  of  elements 

448  7.  ang  -  Look  direction  in  degs:  0  deg  endfire;  90  deg  broadside 

449  7.  180  opposite  endfire 

450  7.  elem  -  Element  spacing  in  wavelenths 

451 

452  if  nargin  <  4; 

453 

454  7.  v  =  exp(  j*2*pi*elem*  [0 :  N—  1  ]  ’  *cos  (ang(  :  )  1  *pi/180) ) ; 

455 

456  Cs  =  -j*2*pi*elem* [0:N-1] ’*sin(ang( : ) ’♦pi/lSO) ; 

457  Cc  =  -j*2*pi*elem*  [0:N-1] ’*cos(ang( : ) ’*pi/180) ; 

458  v  =  [Cc  +  Cs  .  *  Cs]  .  *  exp( j*2*pi*elem* [0 : N— 1] ’*cos(ang( : ) ’*pi/180)) ; 

459  else; 

460  Cs  =  -j*2*pi*elem* ( [0 : N— 1] ’  -  (N-l)/2  ) *sin(ang( : ) ’ *pi/180) ; 

461  Cc  =  -j*2*pi*elem* (  [0:N-1] ’  -  (N-l)/2  )*cos (ang( : ) ’ *pi/180) ; 

462  v  =  [Cc  +  Cs  . *  Cs]  . *  exp(j*2*pi*elem* ( [0:N-1] ’  -  (N-l)/2  ) *cos (ang( : ) y *pi/180) ) ; 

463  end; 

464  return; 

465  7. - 

466  f unct ion [v] =van_p3 (N , ang , elem , f lg) 

467 

468  7.  f  unct  ion  [v]  =van_p3(N ,  ang,  elem,  fig) 

469  7. 

470  7.  This  function  computes  the  third  derivative  of  Vandermonde  steering  vectors. 

471  7. 

472  7.  N  -  Number  of  elements 

473  7.  ang  -  Look  direction  in  degs:  0  deg  endfire;  90  deg  broadside 

474  7.  180  opposite  endfire 

475  7.  elem  -  Element  spacing  in  wavelenths 

476 

477  if  nargin  <  4; 

478  Cs  =  -j*2*pi*elem* [0:N-1] ’♦sinCangC : ) ’♦pi/lSO) ; 

479  Cc  =  -j*2*pi*elem*  [0:N-1] ,*cos(ang( : ) ’*pi/180) ; 

480  eo  =  exp(j*2*pi*elem* [0:N-1] ’*cos(ang( : ) ’*pi/180)) ; 
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481  v  =  (  -Cs  +  Cc.*Cs  +  2*Cs.*Cc  +  Cs.~3  )  . *  eo; 

482  else; 

483  Cs  =  -j*2*pi*elem* ( [0 : N— 1] 1  -  (N-l)/2  ) *sin(ang( : ) ’ *pi/180) ; 

484  Cc  =  -j*2*pi*elem* ( [0 : N— 1] ’  -  (N-l)/2  ) *cos (ang( : ) ’ *pi/180) ; 

485  eo  =  exp( j*2*pi*elem* ( [0 : N-l] *  -  (N-l)/2  ) ♦cos (ang( : ) ’ ♦pi/180) ) ; 

486  v  =  (  -Cs  +  Cc.^Cs  +  2*Cs.*Cc  +  Cs.~3  )  .+60; 

487  end; 

488  return; 

489  7. - 

490  function [R,Ri] =getRvan(N , wl_dB , tht ,thl_dB,d,f lg) ; 

491 

492  7,  function  [R,Ri]  =getRvan(N  , wl.dB ,  tht  ,thl_dB  ,d,  fig) ; 

493  7. 

494  7.  This  function  computes  the  data  covariance  for  simple  mutually 

495  7.  incoherent  planewave  signals  in  white  noise. 

496  7. 

497  7.  N  -  No.  of  elements 

498  7.  wl_dB  -  noise  floor  level  in  dB 

499  7.  tht  -  vector  containing  jammer  directions 

500  7.  thl.dB  -  vector  containing  jammer-to-noise  levels  in  dB 

501  7.  d  -  Element  spacing  in  wavelenths 

502  7.  fig  -  Flag  to  make  array  center  the  phase  reference;  otherwise, 

503  7.  it  *  s  the  f  irst  element 

504 

505  R  =  eye(N)*(  1(T  (wl_dB/10)  )  ; 

506  if  nargin  >  2  &  (  length (tht) ~=0  &  length (thl_dB) ~=0  ); 

507  7»f  unction  [v]  =van  (N ,  ang ,  elem) 

508  if  exist ( ’fig’ ) ;  V  =  van(N,tht ,d,f lg) ;  else;  V  =  van(N , tht ,d) ;  end; 

509  if  length(V) >0 ;  R  =  R  +  V+diag(  10. ~(  (  thl.dB  +  wl.dB)  /  10)  )+V>; 

510  R  =  hm(  R  );  end; 

511  end;  7.  if  nargin  >  2  &  (  length  (tht )  ~=0  &  length (thl.dB) ~=0  ); 

512 

513  if  nargin  >  3  &  (length (tht ) ~=0  &  length (thl.dB) ~=0) ; 

514  Ri  =  10~(  -wl.dB/10  )♦(  eye(N)  -  ... 

515  hm(  V  ♦  hm(  inv(  V’*V+. . . 

516  10~ (wl.dB/ 10) *diag(  10. ~(  -(  thl.dB  +  wl.dB)  /  10)  ) 

517  elseif  nargin  >  3  &  (  length (tht )==0  I  length(thl.dB) ==0  );  ; 

518  Ri  =  10~(  -wl.dB/10  )♦  eye(N)  ; 

519  elseif  nargin  ==  2;  Ri  =  10~(  -wl.dB/10  )*  eye(N)  ; 

520  end; 

521  return; 

522  7. - 

523  function  [Rsqrt ,Rsqrti ,DetRsqrt]  =  get.Rsqrt (R) ; 

524 

525  7»  function  [Rsqrt , Rsqrt i  ,DetRsqrt]  =  get.Rsqrt (R)  ; 

526  7. 

527  7.  This  code  uses  the  eigen-decomposition  to  obtain  a  symmetric  matrix 

528  7#  square  root . 

529  7. 

530  7*  Author:  Christ  D.  Richmond 

531  7.  Date:  11/24/04 

532 

533  [Qr,Ir]  =  eig(R);  Ird  =  diag(Ir) ;  ir.idx  =  f ind(Ird) ;  Iir  =  Ird; 

534  Iir (ir.idx)  =  sqrt(  1 . /Ird(ir.idx)  ); 
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)  )*v>  )  ) 


535  Rsqrt  =  hm(  Qr  *  diag(  sqrt(  Ird(ir_idx)  )  )  *  Qr’  ); 

536  Rsqrti  =  hm(  Qr  *  diag(  Iir  )  *  Qr’  ); 

537  DetRsqrt  =  prod(  Ird  );  %  Determinant 

538  return; 

539  7. - 

540  function [] =vln(vll , varargin) ; 

541 

542  '/,  function [] =vln( vll , varargin) ; 

543  7. 

544  7.  This  plotting  function  plots  vertical  lines  at  locations  indicated  by  vector  vll. 

545  7. 

546  7.  vll  -  vertical  line  locations 

547  7. 

548  7.  Author:  Christ  D.  Richmond 

549  7. 

550  figure(gcf);  ax  =  axis;  hold  on; 

551  for  u  =  1 : length(vll) ;  vo  =  vll(u);  xx  =  vo*[l  1];  yy  =  ax(3:4); 

552  plot (xx,yy , varargin{ : >) ;  end;  hold  off; 

553  7. - 

554  function  [out]  =  MVDR.f nc (param.val) ; 

555 

556  7*  function  [out]  =  MVDR.fnc (param.val) ; 

557  7. 

558  7.  The  MVDR  /  Capon  Spectral  Estimator  parameterized  as  a  function  of 

559  7.  the  unknown  signal  parameters  (angle) 

560  7. 

561  7.  Author:  Christ  D.  Richmond 

562 

563  global  run.params 

564 

565  7.  f unction [v]  =van(N, ang, elem, f  lg)  ; 

566  v_a  =  van (run_params . N ,  param.val , run.params . d) ;  V  =  run.params . Ri.hat  *  v_a; 

567  den  =  (  abs(  sum(  conj (  v_a  )  . *  V,  1  )  )  ); 

568  out  =  den;  7.  Written  upside  down  because  matlab  search  routine  finds  "MINIMA" 

569  7. - 

570  function  [F]  =  cumF(N,M,d,x) 

571 

572  7.  function [F]  =  cumF(N,M,d,x) 

573  7. 

574  7.  This  function  computes  the  cumulative  distribution  function  for  the  complex  central  F. 

575  7. 

576  7.  N  and  M  are  deg  of  freed,  for  complex  non-central  F_{N,M}(d) 

577  7.  d  (  =  delta~2  of  thesis  p.50)  is  non-centrality  parameter 

578  7.  x  value  at  which  to  evaluate  the  cdf  of  F_{N,M}(d) 

579  7. 

580  7.  Author:  Christ  D.  Richmond 

581  7. 

582 

583  idxo  =  find(  x==0  );  7*  idx  =  setdiff  (1 :  length  (x)  ,  idx)  ; 

584  x(idxo)  =  1;  7.  anything  but  zero... 

585 

586  s=0;  for  k=0:M-l;  s  =  s  +  bcf 2CN+M-1 ,k+N) * (x) . ~k . *Gam(d*l . / (1+x) ,k+l) ;  end; 

587  F  =  exp(  N  *  log(x)  -  (N+M-l)  *  log(l  +  x)  +  log(s)  );  F(idxo)  =  0; 

588 
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589  function [g] =Gam(y,m) ;  g  =  1-gammainc (y ,m) ;  return; 

590  function  [y]  =  bcf2(n,k); 

591  for  u=l : length (k) ; 

592  y(u)  =  exp(  gammaln(n+l)  -  gammaln(k(u)+l)  -  gammaln( (n-k(u) )  +  l)  ); 

593  end;  return; 

594  % - 

595  f unct ion [perc] =amidon(t0, iter , total) ; 

596 

597  tnew  =  clock;  startmin  =  60*t0(4)  +  t0(5)  +  t0(6)/60; 

598  pastmin  =  60*tnew(4)  +  tnew (5)  +  tnew(6)/60; 

599  sofarmin  =  pastmin  -  startmin;  min_iter  =  sof armin/iter ; 

600  min_left  =  min_iter  *  (  total  -  iter  ) ; 

601  disp(  [num2str  (100*iter/total)  *'/,  done;  '  num2str  (min.lef  t)  }  Minutes,  or  *  ... 

602  num2str (min.lef t/60)  ’  Hours,  or  ’  num2str (min.lef t/60/24)  ’  Days  left’]); 

603  •/. - 
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ACRONYMS 


AWGN 

Additive  White  Gaussian  Noise 

CBF 

Conventional  Beamforming 

cdf 

Cumulative  Distribution  Function 

CRB 

Cramer-Rao  Bound 

DFT 

Discrete  Fourier  Transform 

DOA 

Direction-of- Arrival 

FFT 

Fast  Fourier  Transform 

FM 

Frequency  Modulation 

HN 

Hawkes-Nehorai  asymptotic  local  error  MSE  prediction  of  Bartlett  algorithm 

hpd 

hermitian  positive  definite 

IE 

Interval  Errors 

KW 

Kaveli  and  Wang 

MFP 

Matched  Field  Processing 

MIE 

Method  of  Interval  Errors 

ML 

Maximum  Likelihood 

MSE 

Mean  Squared  Error 

MUSIC 

Multiple  Signal  Classification 

MVDR 

Minimum  Variance  Distortionless  Response  (adaptive  beamformer) 

NIE 

No  Interval  Errors 

OSF 

Objective  Search  Function 

PCM 

Pulse  Code  Modulation 

pdf 

Probability  Density  Function 

RB 

Rife  and  Boorstyn 

SNR 

Signal-to-Noise  Ratio 

SVD 

Singular  Value  Decomposition 

TKV 

Tufts.  Kot,  and  Vaccaro 
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ACRONYMS  CONTINUED 


UB  Union  Bound 

ULA  Uniform  Linear  Array 

VB  Vaidyanathan-Buckley  asymptotic  local  error  MSE  prediction  of  Capon  algorithm 

WWB  Weiss- Weinstein  Bound 
ZZB  Ziv-Zakai  Bound 
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